feat: separate span power from tx power

gnpy currently uses the same parameter for tx output power and span
input power: this prevents from modelling low tx power effect.
This patch introduces a new tx-cannel-power and uses it to
propagate in ROADM.

Signed-off-by: EstherLerouzic <[email protected]>
Change-Id: Id3ac75e2cb617b513bdb38b51a52e05d15af46f5

gnpy/core/elements.py

@@ -29,7 +29,7 @@
 from logging import getLogger
 
 from gnpy.core.utils import lin2db, db2lin, arrange_frequencies, snr_sum, per_label_average, pretty_summary_print, \
-    watt2dbm, psd2powerdbm
+    watt2dbm, psd2powerdbm, calculate_absolute_min_or_zero
 from gnpy.core.parameters import RoadmParams, FusedParams, FiberParams, PumpParams, EdfaParams, EdfaOperational, \
     RoadmPath, RoadmImpairment
 from gnpy.core.science_utils import NliSolver, RamanSolver
@@ -95,6 +95,7 @@ def __init__(self, *args, **kwargs):
         self.penalties = {}
         self.total_penalty = 0
         self.propagated_labels = [""]
+        self.tx_power = None
 
     def _calc_cd(self, spectral_info):
         """Updates the Transceiver property with the CD of the received channels. CD in ps/nm.
@@ -194,6 +195,7 @@ def __str__(self):
         osnr_ase = per_label_average(self.osnr_ase, self.propagated_labels)
         osnr_ase_01nm = per_label_average(self.osnr_ase_01nm, self.propagated_labels)
         snr_01nm = per_label_average(self.snr_01nm, self.propagated_labels)
+        tx_power_dbm = per_label_average(watt2dbm(self.tx_power), self.propagated_labels)
         cd = mean(self.chromatic_dispersion)
         pmd = mean(self.pmd)
         pdl = mean(self.pdl)
@@ -207,7 +209,8 @@ def __str__(self):
                             f'  CD (ps/nm):                {cd:.2f}',
                             f'  PMD (ps):                  {pmd:.2f}',
                             f'  PDL (dB):                  {pdl:.2f}',
-                            f'  Latency (ms):              {latency:.2f}'])
+                            f'  Latency (ms):              {latency:.2f}',
+                            f'  Actual pch out (dBm):      {pretty_summary_print(tx_power_dbm)}'])
 
         cd_penalty = self.penalties.get('chromatic_dispersion')
         if cd_penalty is not None:
@@ -222,6 +225,7 @@ def __str__(self):
         return result
 
     def __call__(self, spectral_info):
+        self.tx_power = spectral_info.tx_power
         self._calc_snr(spectral_info)
         self._calc_cd(spectral_info)
         self._calc_pmd(spectral_info)
@@ -244,6 +248,8 @@ def __init__(self, *args, params=None, **kwargs):
         # on the path, since it depends on the equalization definition on the degree.
         self.ref_pch_out_dbm = None
         self.loss = 0  # auto-design interest
+        self.loss_pch_db = None
+
         # Optical power of carriers are equalized by the ROADM, so that the experienced loss is not the same for
         # different carriers. The ref_effective_loss records the loss for a reference carrier.
         self.ref_effective_loss = None
@@ -315,11 +321,13 @@ def __str__(self):
             return f'{type(self).__name__} {self.uid}'
 
         total_pch = pretty_summary_print(per_label_average(self.pch_out_dbm, self.propagated_labels))
+        total_loss = pretty_summary_print(per_label_average(self.loss_pch_db, self.propagated_labels))
         return '\n'.join([f'{type(self).__name__} {self.uid}',
-                          f'  type_variety:            {self.type_variety}',
-                          f'  effective loss (dB):     {self.ref_effective_loss:.2f}',
-                          f'  reference pch out (dBm): {self.ref_pch_out_dbm:.2f}',
-                          f'  actual pch out (dBm):    {total_pch}'])
+                          f'  Type_variety:            {self.type_variety}',
+                          f'  Reference loss (dB):     {self.ref_effective_loss:.2f}',
+                          f'  Actual loss (dB):        {total_loss}',
+                          f'  Reference pch out (dBm): {self.ref_pch_out_dbm:.2f}',
+                          f'  Actual pch out (dBm):    {total_pch}'])
 
     def get_roadm_target_power(self, spectral_info: SpectralInformation = None) -> Union[float, ndarray]:
         """Computes the power in dBm for a reference carrier or for a spectral information.
@@ -380,9 +388,13 @@ def propagate(self, spectral_info, degree, from_degree):
         There is no difference for add or express : the same target is applied.
         For the moment propagate operates with spectral info carriers all having the same source or destination.
         """
+        # record input powers to compute the actual loss at the end of the process
+        input_power_dbm = watt2dbm(spectral_info.signal + spectral_info.nli + spectral_info.ase)
         # apply min ROADM loss if it exists
         roadm_maxloss_db = self.get_roadm_path(from_degree, degree).impairment.maxloss
         spectral_info.apply_attenuation_db(roadm_maxloss_db)
+        # records the total power after applying minimum loss
+        net_input_power_dbm = watt2dbm(spectral_info.signal + spectral_info.nli + spectral_info.ase)
         # find the target power for the reference carrier
         ref_per_degree_pch = self.get_per_degree_ref_power(degree)
         # find the target powers for each signal carrier
@@ -392,15 +404,19 @@ def propagate(self, spectral_info, degree, from_degree):
         # Depending on propagation upstream from this ROADM, the input power might be smaller than
         # the target power out configured for this ROADM degree's egress. Since ROADM does not amplify,
         # the power out of the ROADM for the ref channel is the min value between target power and input power.
-        # (TODO add a minimum loss for the ROADM crossing)
-        self.ref_pch_out_dbm = min(self.ref_pch_in_dbm[from_degree] - roadm_maxloss_db, ref_per_degree_pch)
+        ref_pch_in_dbm = self.ref_pch_in_dbm[from_degree]
+        # Calculate the output power for the reference channel (only for visualization)
+        self.ref_pch_out_dbm = min(ref_pch_in_dbm - roadm_maxloss_db, ref_per_degree_pch)
+
         # Definition of effective_loss:
         # Optical power of carriers are equalized by the ROADM, so that the experienced loss is not the same for
         # different carriers. effective_loss records the loss for the reference carrier.
-        self.ref_effective_loss = self.ref_pch_in_dbm[from_degree] - self.ref_pch_out_dbm
-        input_power = spectral_info.signal + spectral_info.nli + spectral_info.ase
+        # Calculate the effective loss for the reference channel
+        self.ref_effective_loss = ref_pch_in_dbm - self.ref_pch_out_dbm
+
+        # Calculate the target power per channel according to the equalization policy
         target_power_per_channel = per_degree_pch + spectral_info.delta_pdb_per_channel
-        # Computation of the per channel target power according to equalization policy
+        # Computation of the correction according to equalization policy
         # If target_power_per_channel has some channels power above input power, then the whole target is reduced.
         # For example, if user specifies delta_pdb_per_channel:
         # freq1: 1dB, freq2: 3dB, freq3: -3dB, and target is -20dBm out of the ROADM,
@@ -415,17 +431,25 @@ def propagate(self, spectral_info, degree, from_degree):
         # that had the min power.
         # This change corresponds to a discussion held during coders call. Please look at this document for
         # a reference: https://telecominfraproject.atlassian.net/wiki/spaces/OOPT/pages/669679645/PSE+Meeting+Minutes
-        correction = (abs(watt2dbm(input_power) - target_power_per_channel)
-                      - (watt2dbm(input_power) - target_power_per_channel)) / 2
+        correction = calculate_absolute_min_or_zero(net_input_power_dbm - target_power_per_channel)
         new_target = target_power_per_channel - correction
-        delta_power = watt2dbm(input_power) - new_target
+        delta_power = net_input_power_dbm - new_target
 
         spectral_info.apply_attenuation_db(delta_power)
-        spectral_info.pmd = sqrt(spectral_info.pmd ** 2
-                                 + self.get_roadm_path(from_degree=from_degree, to_degree=degree).impairment.pmd ** 2)
-        spectral_info.pdl = sqrt(spectral_info.pdl ** 2
-                                 + self.get_roadm_path(from_degree=from_degree, to_degree=degree).impairment.pdl ** 2)
+
+        # Update the PMD information
+        pmd_impairment = self.get_roadm_path(from_degree=from_degree, to_degree=degree).impairment.pmd
+        spectral_info.pmd = sqrt(spectral_info.pmd ** 2 + pmd_impairment ** 2)
+
+        # Update the PMD information
+        pdl_impairment = self.get_roadm_path(from_degree=from_degree, to_degree=degree).impairment.pdl
+        spectral_info.pdl = sqrt(spectral_info.pdl ** 2 + pdl_impairment ** 2)
+
+        # Update the per channel power with the result of propagation
         self.pch_out_dbm = watt2dbm(spectral_info.signal + spectral_info.nli + spectral_info.ase)
+
+        # Update the loss per channel and the labels
+        self.loss_pch_db = input_power_dbm - self.pch_out_dbm
         self.propagated_labels = spectral_info.label
 
     def set_roadm_paths(self, from_degree, to_degree, path_type, impairment_id=None):

gnpy/core/info.py

@@ -52,7 +52,7 @@ class SpectralInformation(object):
 
     def __init__(self, frequency: array, baud_rate: array, slot_width: array, signal: array, nli: array, ase: array,
                  roll_off: array, chromatic_dispersion: array, pmd: array, pdl: array, latency: array,
-                 delta_pdb_per_channel: array, tx_osnr: array, label: array):
+                 delta_pdb_per_channel: array, tx_osnr: array, tx_power: array, label: array):
         indices = argsort(frequency)
         self._frequency = frequency[indices]
         self._df = outer(ones(frequency.shape), frequency) - outer(frequency, ones(frequency.shape))
@@ -80,6 +80,7 @@ def __init__(self, frequency: array, baud_rate: array, slot_width: array, signal
         self._latency = latency[indices]
         self._delta_pdb_per_channel = delta_pdb_per_channel[indices]
         self._tx_osnr = tx_osnr[indices]
+        self._tx_power = tx_power[indices]
         self._label = label[indices]
 
     @property
@@ -188,6 +189,14 @@ def tx_osnr(self):
     def tx_osnr(self, tx_osnr):
         self._tx_osnr = tx_osnr
 
+    @property
+    def tx_power(self):
+        return self._tx_power
+
+    @tx_power.setter
+    def tx_power(self, tx_power):
+        self._tx_power = tx_power
+
     @property
     def channel_number(self):
         return self._channel_number
@@ -232,6 +241,7 @@ def __add__(self, other: SpectralInformation):
                                        delta_pdb_per_channel=append(self.delta_pdb_per_channel,
                                                                     other.delta_pdb_per_channel),
                                        tx_osnr=append(self.tx_osnr, other.tx_osnr),
+                                       tx_power=append(self.tx_power, other.tx_power),
                                        label=append(self.label, other.label))
         except SpectrumError:
             raise SpectrumError('Spectra cannot be summed: channels overlapping.')
@@ -251,6 +261,7 @@ def create_arbitrary_spectral_information(frequency: Union[ndarray, Iterable, fl
                                           signal: Union[float, ndarray, Iterable],
                                           baud_rate: Union[float, ndarray, Iterable],
                                           tx_osnr: Union[float, ndarray, Iterable],
+                                          tx_power: Union[float, ndarray, Iterable] = None,
                                           delta_pdb_per_channel: Union[float, ndarray, Iterable] = 0.,
                                           slot_width: Union[float, ndarray, Iterable] = None,
                                           roll_off: Union[float, ndarray, Iterable] = 0.,
@@ -277,67 +288,71 @@ def create_arbitrary_spectral_information(frequency: Union[ndarray, Iterable, fl
         ase = zeros(number_of_channels)
         delta_pdb_per_channel = full(number_of_channels, delta_pdb_per_channel)
         tx_osnr = full(number_of_channels, tx_osnr)
+        tx_power = full(number_of_channels, tx_power)
         label = full(number_of_channels, label)
         return SpectralInformation(frequency=frequency, slot_width=slot_width,
                                    signal=signal, nli=nli, ase=ase,
                                    baud_rate=baud_rate, roll_off=roll_off,
                                    chromatic_dispersion=chromatic_dispersion,
                                    pmd=pmd, pdl=pdl, latency=latency,
                                    delta_pdb_per_channel=delta_pdb_per_channel,
-                                   tx_osnr=tx_osnr, label=label)
+                                   tx_osnr=tx_osnr, tx_power=tx_power, label=label)
     except ValueError as e:
         if 'could not broadcast' in str(e):
             raise SpectrumError('Dimension mismatch in input fields.')
         else:
             raise
 
 
-def create_input_spectral_information(f_min, f_max, roll_off, baud_rate, power, spacing, tx_osnr, delta_pdb=0):
+def create_input_spectral_information(f_min, f_max, roll_off, baud_rate, spacing, tx_osnr, tx_power,
+                                      delta_pdb=0):
     """Creates a fixed slot width spectral information with flat power.
     all arguments are scalar values"""
     number_of_channels = automatic_nch(f_min, f_max, spacing)
     frequency = [(f_min + spacing * i) for i in range(1, number_of_channels + 1)]
     delta_pdb_per_channel = delta_pdb * ones(number_of_channels)
     label = [f'{baud_rate * 1e-9 :.2f}G' for i in range(number_of_channels)]
-    return create_arbitrary_spectral_information(frequency, slot_width=spacing, signal=power, baud_rate=baud_rate,
+    return create_arbitrary_spectral_information(frequency, slot_width=spacing, signal=tx_power, baud_rate=baud_rate,
                                                  roll_off=roll_off, delta_pdb_per_channel=delta_pdb_per_channel,
-                                                 tx_osnr=tx_osnr, label=label)
+                                                 tx_osnr=tx_osnr, tx_power=tx_power, label=label)
 
 
 def carriers_to_spectral_information(initial_spectrum: dict[float, Carrier],
                                      power: float) -> SpectralInformation:
     """Initial spectrum is a dict with key = carrier frequency, and value a Carrier object.
     :param initial_spectrum: indexed by frequency in Hz, with power offset (delta_pdb), baudrate, slot width,
-    tx_osnr and roll off.
+    tx_osnr, tx_power and roll off.
     :param power: power of the request
     """
     frequency = list(initial_spectrum.keys())
-    signal = [power * db2lin(c.delta_pdb) for c in initial_spectrum.values()]
+    signal = [c.tx_power for c in initial_spectrum.values()]
     roll_off = [c.roll_off for c in initial_spectrum.values()]
     baud_rate = [c.baud_rate for c in initial_spectrum.values()]
     delta_pdb_per_channel = [c.delta_pdb for c in initial_spectrum.values()]
     slot_width = [c.slot_width for c in initial_spectrum.values()]
     tx_osnr = [c.tx_osnr for c in initial_spectrum.values()]
+    tx_power = [c.tx_power for c in initial_spectrum.values()]
     label = [c.label for c in initial_spectrum.values()]
     p_span0 = watt2dbm(power)
     return create_arbitrary_spectral_information(frequency=frequency, signal=signal, baud_rate=baud_rate,
                                                  slot_width=slot_width, roll_off=roll_off,
                                                  delta_pdb_per_channel=delta_pdb_per_channel, tx_osnr=tx_osnr,
-                                                 label=label)
+                                                 tx_power=tx_power, label=label)
 
 
 @dataclass
 class Carrier:
     """One channel in the initial mixed-type spectrum definition, each type being defined by
     its delta_pdb (power offset with respect to reference power), baud rate, slot_width, roll_off
-    and tx_osnr. delta_pdb offset is applied to target power out of Roadm.
+    tx_power, and tx_osnr. delta_pdb offset is applied to target power out of Roadm.
     Label is used to group carriers which belong to the same partition when printing results.
     """
     delta_pdb: float
     baud_rate: float
     slot_width: float
     roll_off: float
     tx_osnr: float
+    tx_power: float
     label: str
 
 

gnpy/core/network.py

@@ -218,7 +218,7 @@ def estimate_raman_gain(node, equipment, power_dbm):
             # do not compute twice to save on time
             return node.estimated_gain
         spectral_info = create_input_spectral_information(f_min=f_min, f_max=f_max, roll_off=roll_off,
-                                                          baud_rate=baud_rate, power=power, spacing=spacing,
+                                                          baud_rate=baud_rate, tx_power=power, spacing=spacing,
                                                           tx_osnr=tx_osnr)
         pin = watt2dbm(sum(spectral_info.signal))
         attenuation_in_db = node.params.con_in + node.params.att_in
@@ -304,9 +304,11 @@ def set_egress_amplifier(network, this_node, equipment, pref_ch_db, pref_total_d
         prev_node = this_node
         node = oms
         if isinstance(this_node, elements.Transceiver):
-            # for the time being use the same power for the target of roadms and for transceivers
-            # TODO: This should be changed when introducing a power parameter dedicated to transceivers
-            this_node_out_power = pref_ch_db
+            # todo change pref to a ref channel
+            if equipment['SI']['default'].tx_power_dbm is not None:
+                this_node_out_power = equipment['SI']['default'].tx_power_dbm
+            else:
+                this_node_out_power = pref_ch_db
         if isinstance(this_node, elements.Roadm):
             # get target power out from ROADM for the reference carrier based on equalization settings
             this_node_out_power = this_node.get_per_degree_ref_power(degree=node.uid)

gnpy/core/utils.py

@@ -9,7 +9,7 @@
 """
 
 from csv import writer
-from numpy import pi, cos, sqrt, log10, linspace, zeros, shape, where, logical_and, mean
+from numpy import pi, cos, sqrt, log10, linspace, zeros, shape, where, logical_and, mean, array
 from scipy import constants
 from copy import deepcopy
 
@@ -452,3 +452,20 @@ def restore_order(elements, order):
     [3, 2, 7]
     """
     return [elements[i[0]] for i in sorted(enumerate(order), key=lambda x:x[1]) if elements[i[0]] is not None]
+
+
+def calculate_absolute_min_or_zero(x: array) -> array:
+    """Calculates the element-wise absolute minimum between the x and zero.
+
+    Parameters:
+    x (array): The first input array.
+
+    Returns:
+    array: The element-wise absolute minimum between x and zero.
+
+    Example:
+    >>> x = array([-1, 2, -3])
+    >>> calculate_absolute_min_or_zero(x)
+    array([1., 0., 3.])
+    """
+    return (abs(x) - x) / 2

gnpy/example-data/eqpt_config.json

@@ -290,6 +290,7 @@
             "spacing": 50e9,
             "power_dbm": 0,
             "power_range_db": [0, 0, 1],
+            "tx_power_dbm": 0,
             "roll_off": 0.15,
             "tx_osnr": 40,
             "sys_margins": 2

gnpy/tools/cli_examples.py

@@ -194,8 +194,10 @@ def transmission_main_example(args=None):
     params['effective_freq_slot'] = None
     trx_params = trx_mode_params(equipment)
     trx_params['power'] = dbm2watt(equipment['SI']['default'].power_dbm)
+    trx_params['tx_power'] = dbm2watt(equipment['SI']['default'].power_dbm)
     if args.power:
         trx_params['power'] = dbm2watt(float(args.power))
+        trx_params['tx_power'] = dbm2watt(float(args.power))
     params.update(trx_params)
     initial_spectrum = None
     params['nb_channel'] = automatic_nch(trx_params['f_min'], trx_params['f_max'], trx_params['spacing'])
@@ -372,7 +374,8 @@ def path_requests_run(args=None):
         'effective_freq_slot': None,
         'nb_channel': automatic_nch(equipment['SI']['default'].f_min, equipment['SI']['default'].f_max,
                                     equipment['SI']['default'].spacing),
-        'power': dbm2watt(equipment['SI']['default'].power_dbm)
+        'power': dbm2watt(equipment['SI']['default'].power_dbm),
+        'tx_power': dbm2watt(equipment['SI']['default'].power_dbm)
     }
     trx_params = trx_mode_params(equipment)
     params.update(trx_params)