Optional Flows in LinearConverter

[baumbude]

optional and obligatory flows

[FBumann]

Staring a more generic approach:

Making flows optional would make modeling certain things much easier.
Choosing the right input flow or changing between different possible inputs while having the same output would be possible.
Further, things like converting to different temperature levels could also be modeled.

[FBumann]

For a nice user experience, it would be nice to have such a possible workflow:

import numpy as np
import fx

# Define inputs
ins = [fx.Flow(label='El Input', bus=el_bus)]

# Define outputs
outs = [
    fx.Flow(label='Thermal Output 1', bus=heat_low_temp, size=100),
    fx.Flow(label='Thermal Output 2', bus= heat_mid_temp, size=100),
    fx.Flow(label='Thermal Output 3', bus= heat_high_temp, size=100)
]

# Define multiple conversion factors. One of them is active at a time. Which one is part of the optimization
convs = [
    {ins[0]: 3.5, outs[0]: 1},
    {ins[0]: 3.0, outs[1]: 1},
    {ins[0]: 2.5, outs[2]: 1}
]

# Define the Heat Pump
boiler = fx.LinearConverter('Heat Pump', inputs=ins, outputs=outs, conversion_factors=convs)

[FBumann]

Mathematically:

Controlling the Number of Active Equations in a LinearConverterModel

1. Binary Variables for each conversion factor: $z_n$ (1 if the $n$-th equation is active, 0 otherwise).
2. Equations (Big-M method):
   • Change the Equation for each conversion factor $n$ to an Inequation, using a binary variable to relax the equation

$$ \sum_{I}a_i x_i - \sum_{J} b_j y_j \leq M \cdot (1 - z_n) $$

• When $z_i = 1$, equation $a_i x = b_i y$ is enforced.
• When $z_i = 0$, equation is relaxed.

4. Control the Number of Active Equations at a time:
   • Exactly 1 active: $\sum_n z_n = 1$
5. **Make sure all flow rates $v$ not defined by the active conversion_factor are 0
   • $v \leq M \cdot (1-z_n)$

Issues with other functionality:

• Flows that can not be 0 would lead to not solvable models

This was tested with a super basic example and worked

# -*- coding: utf-8 -*-
"""
Energy System Optimization Example
Developed by Felix Panitz and Peter Stange
Chair of Building Energy Systems and Heat Supply, Technische Universität Dresden
"""

import numpy as np
import flixOpt as fx

if __name__ == '__main__':

    # --- Define Thermal Load Profile ---
    # Load profile (e.g., kW) for heating demand over time
    thermal_load_profile = np.array([30., 0., 20.])
    datetime_series = fx.create_datetime_array('2020-01-01', 3, 'h')

    # --- Define Energy Buses ---
    # These represent the different energy carriers in the system
    electricity_bus = fx.Bus('Electricity')
    heat_bus = fx.Bus('District Heating')
    fuel_bus = fx.Bus('Natural Gas')

    # --- Define Objective Effect (Cost) ---
    # Cost effect representing the optimization objective (minimizing costs)
    cost_effect = fx.Effect('costs', '€', 'Cost', is_standard=True, is_objective=True)

    # --- Define Flow System Components ---
    # Boiler component with thermal output (heat) and fuel input (gas)
    ins = [fx.Flow(label='Fuel Input', bus=fuel_bus)]
    outs= [fx.Flow(label=f'Thermal Output {i}', bus=heat_bus, size=100,
                   relative_maximum= np.array([1 if j == i-1 else 0 for j in range(3)]),
                   effects_per_flow_hour=1) for i in range(1,4)]
    convs = [{ins[0]: 0.8, flow: 1} for flow in outs]
    boiler = fx.LinearConverter('Boiler', inputs=ins, outputs=outs,
                                conversion_factors=convs)

    # Heat load component with a fixed thermal demand profile
    heat_load = fx.Sink('Heat Demand',
        sink=fx.Flow(label='Thermal Load', bus=heat_bus, size=1, fixed_relative_profile=thermal_load_profile))

    # Gas source component with cost-effect per flow hour
    gas_source = fx.Source('Natural Gas Tariff',
        source=fx.Flow(label='Gas Flow', bus=fuel_bus, size=1000, effects_per_flow_hour=0.04))  # 0.04 €/kWh

    # --- Build the Flow System ---
    # Add all components and effects to the system
    flow_system = fx.FlowSystem(datetime_series)
    flow_system.add_elements(cost_effect, boiler, heat_load, gas_source)

    # --- Define and Run Calculation ---
    calculation = fx.FullCalculation('Simulation1', flow_system)
    calculation.do_modeling()

    # --- Solve the Calculation and Save Results ---
    calculation.solve(fx.solvers.HighsSolver(), save_results=True)

    # --- Load and Analyze Results ---
    # Load results and plot the operation of the District Heating Bus
    results = fx.results.CalculationResults(calculation.name, folder='results')
    results.plot_operation('District Heating', 'area')

    # Print results to the console. Check Results in file or perform more plotting
    from pprint import pprint
    pprint(calculation.results())
    print(f'Look into .yaml and .json file for results')

Results:

{'Buses': {'District Heating': {'excess_input': array([0., 0., 0.]),
                                'excess_output': array([0., 0., 0.])},
           'Natural Gas': {'excess_input': array([0., 0., 0.]),
                           'excess_output': array([ 0., -0.,  0.])}},
 'Components': {'Boiler': {'Fuel Input': {'flow_rate': array([37.5,  0. , 25. ]),
                                          'sumFlowHours': 62.5},
                           'Thermal Output 1': {'flow_rate': array([30.,  0.,  0.]),
                                                'sumFlowHours': 30.0},
                           'Thermal Output 2': {'flow_rate': array([0., 0., 0.]),
                                                'sumFlowHours': 0.0},
                           'Thermal Output 3': {'flow_rate': array([ 0.,  0., 20.]),
                                                'sumFlowHours': 20.0},
                           'conversion_factor_0': array([1, 1, 0], dtype=int8),
                           'conversion_factor_1': array([0, 0, 0], dtype=int8),
                           'conversion_factor_2': array([0, 0, 1], dtype=int8)},
                'Heat Demand': {'Thermal Load': {'flow_rate': array([30.,  0., 20.]),
                                                 'sumFlowHours': 50.0}},
                'Natural Gas Tariff': {'Gas Flow': {'flow_rate': array([37.5,  0. , 25. ]),
                                                    'sumFlowHours': 62.5}}},
 'Effects': {'costs': {'all': {'Shares': {'invest': -0.0, 'operation': 52.5},
                               'all_sum': 52.5},
                       'invest': {'Shares': {}, 'invest_sum': -0.0},
                       'operation': {'Shares': {'Boiler__Thermal Output 1__effects_per_flow_hour': array([30.,  0.,  0.]),
                                                'Boiler__Thermal Output 2__effects_per_flow_hour': array([0., 0., 0.]),
                                                'Boiler__Thermal Output 3__effects_per_flow_hour': array([ 0.,  0., 20.]),
                                                'Natural Gas Tariff__Gas Flow__effects_per_flow_hour': array([1.5, 0. , 1. ])},
                                     'operation_sum': 52.5,
                                     'operation_sum_TS': array([31.5,  0. , 21. ])}},
             'penalty': {'Shares': {'District Heating__excess_input': 0.0,
                                    'District Heating__excess_output': 0.0,
                                    'Natural Gas__excess_input': 0.0,
                                    'Natural Gas__excess_output': 0.0},
                         'penalty_sum': 0.0}},
 'Objective': 52.5,
 'Others': {},
 'Time': array(['2020-01-01T00', '2020-01-01T01', '2020-01-01T02', '2020-01-01T03'],
      dtype='datetime64[h]'),
 'Time intervals in hours': array([1., 1., 1.])}

  Constraints:
    Components:
      Boiler:
        _self:
        - 'EQ Boiler_use_rates      [1/3]:      1.0 = 1.0 * Boiler_conversion_factor_0[0] + 1.0 * Boiler_conversion_factor_1[0] + 1.0 * Boiler_conversion_factor_2[0]'
        - 'INEQ Boiler_conversion_factor_0_ub [1/3]: 1000000.0 >= 0.8 * Boiler__Fuel Input_flow_rate[0] + -1.0 * Boiler__Thermal Output 1_flow_rate[0] + 1e+06 * Boiler_conversion_factor_0[0]'
        - 'INEQ Boiler_conversion_factor_0_lb [1/3]: 1000000.0 >= -0.8 * Boiler__Fuel Input_flow_rate[0] + 1.0 * Boiler__Thermal Output 1_flow_rate[0] + 1e+06 * Boiler_conversion_factor_0[0]'
        - 'INEQ Boiler_conversion_factor_0_restrict_Thermal Output 3 [1/3]: 1000000.0 >= 1.0 * Boiler__Thermal Output 3_flow_rate[0] + 1e+06 * Boiler_conversion_factor_0[0]'
        - 'INEQ Boiler_conversion_factor_0_restrict_Thermal Output 2 [1/3]: 1000000.0 >= 1.0 * Boiler__Thermal Output 2_flow_rate[0] + 1e+06 * Boiler_conversion_factor_0[0]'
        - 'INEQ Boiler_conversion_factor_1_ub [1/3]: 1000000.0 >= 0.8 * Boiler__Fuel Input_flow_rate[0] + -1.0 * Boiler__Thermal Output 2_flow_rate[0] + 1e+06 * Boiler_conversion_factor_1[0]'
        - 'INEQ Boiler_conversion_factor_1_lb [1/3]: 1000000.0 >= -0.8 * Boiler__Fuel Input_flow_rate[0] + 1.0 * Boiler__Thermal Output 2_flow_rate[0] + 1e+06 * Boiler_conversion_factor_1[0]'
        - 'INEQ Boiler_conversion_factor_1_restrict_Thermal Output 3 [1/3]: 1000000.0 >= 1.0 * Boiler__Thermal Output 3_flow_rate[0] + 1e+06 * Boiler_conversion_factor_1[0]'
        - 'INEQ Boiler_conversion_factor_1_restrict_Thermal Output 1 [1/3]: 1000000.0 >= 1.0 * Boiler__Thermal Output 1_flow_rate[0] + 1e+06 * Boiler_conversion_factor_1[0]'
        - 'INEQ Boiler_conversion_factor_2_ub [1/3]: 1000000.0 >= 0.8 * Boiler__Fuel Input_flow_rate[0] + -1.0 * Boiler__Thermal Output 3_flow_rate[0] + 1e+06 * Boiler_conversion_factor_2[0]'
        - 'INEQ Boiler_conversion_factor_2_lb [1/3]: 1000000.0 >= -0.8 * Boiler__Fuel Input_flow_rate[0] + 1.0 * Boiler__Thermal Output 3_flow_rate[0] + 1e+06 * Boiler_conversion_factor_2[0]'
        - 'INEQ Boiler_conversion_factor_2_restrict_Thermal Output 1 [1/3]: 1000000.0 >= 1.0 * Boiler__Thermal Output 1_flow_rate[0] + 1e+06 * Boiler_conversion_factor_2[0]'
        - 'INEQ Boiler_conversion_factor_2_restrict_Thermal Output 2 [1/3]: 1000000.0 >= 1.0 * Boiler__Thermal Output 2_flow_rate[0] + 1e+06 * Boiler_conversion_factor_2[0]'
        Fuel Input:
          _self:
          - 'EQ Boiler__Fuel Input_sumFlowHours [1/1]:      0.0 = ∑(1.0 * Boiler__Fuel Input_flow_rate[0]+..) + -1.0 * Boiler__Fuel Input_sumFlowHours[0]'
        Thermal Output 1:
          _self:
          - 'EQ Boiler__Thermal Output 1_sumFlowHours [1/1]:      0.0 = ∑(1.0 * Boiler__Thermal Output 1_flow_rate[0]+..) + -1.0 * Boiler__Thermal Output 1_sumFlowHours[0]'
        Thermal Output 2:
          _self:
          - 'EQ Boiler__Thermal Output 2_sumFlowHours [1/1]:      0.0 = ∑(1.0 * Boiler__Thermal Output 2_flow_rate[0]+..) + -1.0 * Boiler__Thermal Output 2_sumFlowHours[0]'
        Thermal Output 3:
          _self:
          - 'EQ Boiler__Thermal Output 3_sumFlowHours [1/1]:      0.0 = ∑(1.0 * Boiler__Thermal Output 3_flow_rate[0]+..) + -1.0 * Boiler__Thermal Output 3_sumFlowHours[0]'

[FBumann]

In contrary to prior discussions, after some testing, "optional flows" are already possible.
They can be defined like this:

ins = [fx.Flow(label='Fuel Input', bus=fuel_bus)]
outs = [fx.Flow(label='Thermal Output', bus=heat_bus, size=50, relative_maximum=np.array([1,1,0])),
        fx.Flow(label='Thermal Output 2', bus=heat_bus, size=50, relative_maximum=np.array([0,1,1]))]
boiler = fx.LinearConverter('Boiler', inputs=ins, outputs=outs,
                            conversion_factors=[{ins[0]: 1, outs[0]: 1/0.8, outs[1]: 1/0.7}])

The resulting math would be:

• $1 \cdot p_\text{fuel} - 1/0.8 \cdot p_\text{th1} - 1/0.7 \cdot p_\text{th2} = 0$

This effectively allows the optimizer to choose between the different thermal outputs.
But is this also possible with more complex relations, with multiple flows and conversion rates?

Full Script for testing:

# -*- coding: utf-8 -*-
"""
Energy System Optimization Example
Developed by Felix Panitz and Peter Stange
Chair of Building Energy Systems and Heat Supply, Technische Universität Dresden
"""

import numpy as np
import flixOpt as fx

if __name__ == '__main__':

    # --- Define Thermal Load Profile ---
    # Load profile (e.g., kW) for heating demand over time
    thermal_load_profile = np.array([30., 0., 20.])
    datetime_series = fx.create_datetime_array('2020-01-01', 3, 'h')

    # --- Define Energy Buses ---
    # These represent the different energy carriers in the system
    electricity_bus = fx.Bus('Electricity')
    heat_bus = fx.Bus('District Heating')
    fuel_bus = fx.Bus('Natural Gas')

    # --- Define Objective Effect (Cost) ---
    # Cost effect representing the optimization objective (minimizing costs)
    cost_effect = fx.Effect('costs', '€', 'Cost', is_standard=True, is_objective=True)

    # --- Define Flow System Components ---
    # Boiler component with thermal output (heat) and fuel input (gas)
    ins = [fx.Flow(label='Fuel Input', bus=fuel_bus)]
    outs = [fx.Flow(label='Thermal Output', bus=heat_bus, size=50, relative_maximum=np.array([1,1,0])),
            fx.Flow(label='Thermal Output 2', bus=heat_bus, size=50, relative_maximum=np.array([0,1,1]))]
    boiler = fx.LinearConverter('Boiler', inputs=ins, outputs=outs,
                                conversion_factors=[{ins[0]: 1, outs[0]: 1/0.8, outs[1]: 1/0.7}])

    # Heat load component with a fixed thermal demand profile
    heat_load = fx.Sink('Heat Demand',
        sink=fx.Flow(label='Thermal Load', bus=heat_bus, size=1, fixed_relative_profile=thermal_load_profile))

    # Gas source component with cost-effect per flow hour
    gas_source = fx.Source('Natural Gas Tariff',
        source=fx.Flow(label='Gas Flow', bus=fuel_bus, size=1000, effects_per_flow_hour=0.04))  # 0.04 €/kWh

    # --- Build the Flow System ---
    # Add all components and effects to the system
    flow_system = fx.FlowSystem(datetime_series)
    flow_system.add_elements(cost_effect, boiler, heat_load, gas_source)

    # --- Define and Run Calculation ---
    calculation = fx.FullCalculation('Simulation1', flow_system)
    calculation.do_modeling()

    # --- Solve the Calculation and Save Results ---
    calculation.solve(fx.solvers.HighsSolver(), save_results=True)

    # --- Load and Analyze Results ---
    # Load results and plot the operation of the District Heating Bus
    results = fx.results.CalculationResults(calculation.name, folder='results')
    results.plot_operation('District Heating', 'area')

    # Print results to the console. Check Results in file or perform more plotting
    from pprintpp import pprint
    pprint(calculation.results())
    print(f'Look into .yaml and .json file for results')

[baumbude]

Ok, da sollten wir überlegen wie wir diese 2 unterschiedlichen Funktionalitäten zusammenführen können.
Ein Grund-Prinzip von flixopt ist: möglichst viel Funktionalität bereits ohne Binärvariablen. D.h. man müsste die Nichtgleichzeitigkeit oben drauf setzen. Dafür bietet sich an das Feature "avoidFlowAtOnce" (alte Bezeichnung) an.
Habe gerade Idee für spannende Funktionserweiterung derselben...
Cool wäre noch wenn wir Segments gleichartig mitdenken...

[FBumann]

Ich denke auch dass wir die Restriktion ENTWEDER/ODER als Zusatz lassen.