# Mixed Integer Optimization: Pumps and Orifices

Note

This example focuses on how to incorporate mixed integer components into a hydraulic model, and assumes basic exposure to RTC-Tools. To start with basics, see Filling a Reservoir.

Note

By default, if you define any integer or boolean variables in the model, RTC-Tools
will switch from IPOPT to BONMIN. You can modify solver options by overriding
the
`solver_options()`

method. Refer to CasADi’s nlpsol interface for a list of supported solvers.

## The Model

For this example, the model represents a typical setup for the dewatering of lowland areas. Water is routed from the hinterland (modeled as discharge boundary condition, right side) through a canal (modeled as storage element) towards the sea (modeled as water level boundary condition on the left side). Keeping the lowland area dry requires that enough water is discharged to the sea. If the sea water level is lower than the water level in the canal, the water can be discharged to the sea via gradient flow through the orifice (or a weir). If the sea water level is higher than in the canal, water must be pumped.

To discharge water via gradient flow is free, while pumping costs money. The control task is to keep the water level in the canal below a given flood warning level at minimum costs. The expected result is that the model computes a control pattern that makes use of gradient flow whenever possible and activates the pump only when necessary.

The model can be viewed and edited using the OpenModelica Connection Editor
program. First load the Deltares library into OpenModelica Connection Editor,
and then load the example model, located at
`<examples directory>\mixed_integer\model\Example.mo`

. The model `Example.mo`

represents a simple water system with the following elements:

a canal segment, modeled as storage element

`Deltares.ChannelFlow.Hydraulic.Storage.Linear`

,a discharge boundary condition

`Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Discharge`

,a water level boundary condition

`Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Level`

,a pump

`Deltares.ChannelFlow.Hydraulic.Structures.Pump`

an orifice modeled as a pump

`Deltares.ChannelFlow.Hydraulic.Structures.Pump`

In text mode, the Modelica model looks as follows (with annotation statements removed):

```
1model Example
2 // Declare Model Elements
3 Deltares.ChannelFlow.Hydraulic.Storage.Linear storage(A=1.0e6, H_b=0.0, HQ.H(min=0.0, max=0.5));
4 Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Discharge discharge;
5 Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Level level;
6 Deltares.ChannelFlow.Hydraulic.Structures.Pump pump;
7 Deltares.ChannelFlow.Hydraulic.Structures.Pump orifice;
8
9 // Define Input/Output Variables and set them equal to model variables
10 input Modelica.SIunits.VolumeFlowRate Q_pump(fixed=false, min=0.0, max=7.0) = pump.Q;
11 input Boolean is_downhill;
12 input Modelica.SIunits.VolumeFlowRate Q_in(fixed=true) = discharge.Q;
13 input Modelica.SIunits.Position H_sea(fixed=true) = level.H;
14 input Modelica.SIunits.VolumeFlowRate Q_orifice(fixed=false, min=0.0, max=10.0) = orifice.Q;
15 output Modelica.SIunits.Position storage_level = storage.HQ.H;
16 output Modelica.SIunits.Position sea_level = level.H;
17equation
18 // Connect Model Elements
19 connect(orifice.HQDown, level.HQ);
20 connect(storage.HQ, orifice.HQUp);
21 connect(storage.HQ, pump.HQUp);
22 connect(discharge.HQ, storage.HQ);
23 connect(pump.HQDown, level.HQ);
24end Example;
```

The five water system elements (storage, discharge boundary condition, water
level boundary condition, pump, and orifice) appear under the `model Example`

statement. The `equation`

part connects these five elements with the help of
connections. Note that `Pump`

extends the partial model `HQTwoPort`

which
inherits from the connector `HQPort`

. With `HQTwoPort`

, `Pump`

can be
connected on two sides. `level`

represents a model boundary condition (model
is meant in a hydraulic sense here), so it can be connected to one other
element only. It extends the `HQOnePort`

which again inherits from the
connector `HQPort`

.

In addition to elements, the input variables `Q_in`

, `H_sea`

, `Q_pump`

,
and `Q_orifice`

are also defined. Because we want to view the water levels in
the storage element in the output file, we also define output
variables `storage_level`

and `sea_level`

. It is usually easiest to set input
and output variables equal to their corresponding model variable in the same line.

To maintain the linearity of the model, we input the Boolean `is_downhill`

as
a way to keep track of whether water can flow by gravity to the sea. This
variable is not used directly in the hydraulics, but we use it later in the
constraints in the python file.

## The Optimization Problem

The python script consists of the following blocks:

Import of packages

Definition of the optimization problem class

Constructor

Objective function

Definition of constraints

Additional configuration of the solver

A run statement

### Importing Packages

For this example, the import block is as follows:

```
1import numpy as np
2from rtctools.optimization.collocated_integrated_optimization_problem import (
3 CollocatedIntegratedOptimizationProblem,
4)
5from rtctools.optimization.csv_mixin import CSVMixin
6from rtctools.optimization.modelica_mixin import ModelicaMixin
```

Note that we are also importing `inf`

from `numpy`

. We will use this later
in the constraints.

### Optimization Problem

Next, we construct the class by declaring it and inheriting the desired parent classes.

```
10class Example(CSVMixin, ModelicaMixin, CollocatedIntegratedOptimizationProblem):
```

Now we define an objective function. This is a class method that returns the value that needs to be minimized. Here we specify that we want to minimize the volume pumped:

```
18 def objective(self, ensemble_member):
19 # Minimize water pumped. The total water pumped is the integral of the
20 # water pumped from the starting time until the stoping time. In
21 # practice, self.integral() is a summation of all the discrete states.
22 return self.integral("Q_pump", ensemble_member=ensemble_member)
```

Constraints can be declared by declaring the `path_constraints()`

method.
Path constraints are constraints that are applied every timestep. To set a
constraint at an individual timestep, define it inside the `constraints`

method.

The orifice `BooleanSubmergedOrifice`

requires special constraints to be set
in order to work. They are implemented below in the `path_constraints()`

method. their parent classes also declare this method, so we call the
`super()`

method so that we don’t overwrite their behaviour.

```
26 def path_constraints(self, ensemble_member):
27 # Call super to get default constraints
28 constraints = super().path_constraints(ensemble_member)
29 M = 2 # The so-called "big-M"
30
31 # Release through orifice downhill only. This constraint enforces the
32 # fact that water only flows downhill.
33 constraints.append(
34 (self.state("Q_orifice") + (1 - self.state("is_downhill")) * 10, 0.0, 10.0)
35 )
36
37 # Make sure is_downhill is true only when the sea is lower than the
38 # water level in the storage.
39 constraints.append(
40 (
41 self.state("H_sea")
42 - self.state("storage.HQ.H")
43 - (1 - self.state("is_downhill")) * M,
44 -np.inf,
45 0.0,
46 )
47 )
48 constraints.append(
49 (
50 self.state("H_sea") - self.state("storage.HQ.H") + self.state("is_downhill") * M,
51 0.0,
52 np.inf,
53 )
54 )
55
56 # Orifice flow constraint. Uses the equation:
57 # Q(HUp, HDown, d) = width * C * d * (2 * g * (HUp - HDown)) ^ 0.5
58 # Note that this equation is only valid for orifices that are submerged
59 # units: description:
60 w = 3.0 # m width of orifice
61 d = 0.8 # m hight of orifice
62 C = 1.0 # none orifice constant
63 g = 9.8 # m/s^2 gravitational acceleration
64 constraints.append(
65 (
66 ((self.state("Q_orifice") / (w * C * d)) ** 2) / (2 * g)
67 + self.state("orifice.HQDown.H")
68 - self.state("orifice.HQUp.H")
69 - M * (1 - self.state("is_downhill")),
70 -np.inf,
71 0.0,
72 )
73 )
74
75 return constraints
```

Finally, we want to apply some additional configuration, reducing the amount of information the solver outputs:

```
78 def solver_options(self):
79 options = super().solver_options()
80 # Restrict solver output
81 solver = options["solver"]
82 options[solver]["print_level"] = 1
83 return options
```

### Run the Optimization Problem

To make our script run, at the bottom of our file we just have to call
the `run_optimization_problem()`

method we imported on the optimization
problem class we just created.

```
87run_optimization_problem(Example)
```

### The Whole Script

All together, the whole example script is as follows:

```
1import numpy as np
2from rtctools.optimization.collocated_integrated_optimization_problem import (
3 CollocatedIntegratedOptimizationProblem,
4)
5from rtctools.optimization.csv_mixin import CSVMixin
6from rtctools.optimization.modelica_mixin import ModelicaMixin
7from rtctools.util import run_optimization_problem
8
9
10class Example(CSVMixin, ModelicaMixin, CollocatedIntegratedOptimizationProblem):
11 """
12 This class is the optimization problem for the Example. Within this class,
13 the objective, constraints and other options are defined.
14 """
15
16 # This is a method that returns an expression for the objective function.
17 # RTC-Tools always minimizes the objective.
18 def objective(self, ensemble_member):
19 # Minimize water pumped. The total water pumped is the integral of the
20 # water pumped from the starting time until the stoping time. In
21 # practice, self.integral() is a summation of all the discrete states.
22 return self.integral("Q_pump", ensemble_member=ensemble_member)
23
24 # A path constraint is a constraint where the values in the constraint are a
25 # Timeseries rather than a single number.
26 def path_constraints(self, ensemble_member):
27 # Call super to get default constraints
28 constraints = super().path_constraints(ensemble_member)
29 M = 2 # The so-called "big-M"
30
31 # Release through orifice downhill only. This constraint enforces the
32 # fact that water only flows downhill.
33 constraints.append(
34 (self.state("Q_orifice") + (1 - self.state("is_downhill")) * 10, 0.0, 10.0)
35 )
36
37 # Make sure is_downhill is true only when the sea is lower than the
38 # water level in the storage.
39 constraints.append(
40 (
41 self.state("H_sea")
42 - self.state("storage.HQ.H")
43 - (1 - self.state("is_downhill")) * M,
44 -np.inf,
45 0.0,
46 )
47 )
48 constraints.append(
49 (
50 self.state("H_sea") - self.state("storage.HQ.H") + self.state("is_downhill") * M,
51 0.0,
52 np.inf,
53 )
54 )
55
56 # Orifice flow constraint. Uses the equation:
57 # Q(HUp, HDown, d) = width * C * d * (2 * g * (HUp - HDown)) ^ 0.5
58 # Note that this equation is only valid for orifices that are submerged
59 # units: description:
60 w = 3.0 # m width of orifice
61 d = 0.8 # m hight of orifice
62 C = 1.0 # none orifice constant
63 g = 9.8 # m/s^2 gravitational acceleration
64 constraints.append(
65 (
66 ((self.state("Q_orifice") / (w * C * d)) ** 2) / (2 * g)
67 + self.state("orifice.HQDown.H")
68 - self.state("orifice.HQUp.H")
69 - M * (1 - self.state("is_downhill")),
70 -np.inf,
71 0.0,
72 )
73 )
74
75 return constraints
76
77 # Any solver options can be set here
78 def solver_options(self):
79 options = super().solver_options()
80 # Restrict solver output
81 solver = options["solver"]
82 options[solver]["print_level"] = 1
83 return options
84
85
86# Run
87run_optimization_problem(Example)
```

## Running the Optimization Problem

Note

An explaination of bonmin behaviour and output goes here.

## Extracting Results

The results from the run are found in `output/timeseries_export.csv`

. Any
CSV-reading software can import it, but this is how results can be plotted using
the python library matplotlib:

```
from datetime import datetime
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
# Import Data
data_path = "../../../examples/mixed_integer/reference_output/timeseries_export.csv"
results = np.recfromcsv(data_path, encoding=None)
# Get times as datetime objects
times = [datetime.strptime(x, "%Y-%m-%d %H:%M:%S") for x in results["time"]]
# Generate Plot
fig, axarr = plt.subplots(2, sharex=True)
axarr[0].set_title("Water Level and Discharge")
# Upper subplot
axarr[0].set_ylabel("Water Level [m]")
axarr[0].plot(times, results["storage_level"], label="Storage", linewidth=2, color="b")
axarr[0].plot(times, results["sea_level"], label="Sea", linewidth=2, color="m")
axarr[0].plot(
times,
0.5 * np.ones_like(times),
label="Storage Max",
linewidth=2,
color="r",
linestyle="--",
)
# Lower Subplot
axarr[1].set_ylabel("Flow Rate [m³/s]")
# add dots to clarify where the decision variables are defined:
axarr[1].scatter(times, results["q_orifice"], linewidth=1, color="g")
axarr[1].scatter(times, results["q_pump"], linewidth=1, color="r")
# add horizontal lines to the left of these dots, to indicate that the value is attained over an
# entire timestep:
axarr[1].step(times, results["q_orifice"], linewidth=2, where="pre", label="Orifice", color="g")
axarr[1].step(times, results["q_pump"], linewidth=1, where="pre", label="Pump", color="r")
# Format bottom axis label
axarr[-1].xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
# Shrink margins
fig.tight_layout()
# Shrink each axis and put a legend to the right of the axis
for i in range(len(axarr)):
box = axarr[i].get_position()
axarr[i].set_position([box.x0, box.y0, box.width * 0.8, box.height])
axarr[i].legend(loc="center left", bbox_to_anchor=(1, 0.5), frameon=False)
plt.autoscale(enable=True, axis="x", tight=True)
# Output Plot
plt.show()
```

(`Source code`

, `svg`

, `png`

)

## Observations

Note that in the results plotted above, the pump runs with a constantly varying throughput. To smooth out the flow through the pump, consider using goal programming to apply a path goal minimizing the derivative of the pump at each timestep. For an example, see the third goal in Declaring Goals.