Controlling a Cascade of Weirs: Local PID Control vs. Optimization

../../_images/lavant-2976738_1280.jpg

Note

This is a more advanced example that implements multi-objective optimization in RTC-Tools. It also capitalizes on the homotopy techniques available in RTC-Tools. If you are a first-time user of RTC-Tools, see Filling a Reservoir.

One of the advantages of using RTC-Tools for control is that it is capable of making decisions that are optimal for the whole network and for all future timesteps within the model time horizon. This is in contrast to local control algorithms such as PID controllers, where the control decision must be made on past states and local information alone. Furthermore, unlike a PID-style controller, RTC-Tools does not have gain parameters that need to be tuned.

This example models a cascading channel system, and compares a local control scheme using PID controllers with the RTC-Tools approach that uses Goal Programming.

The Model

For this example, water is flowing through a multilevel channel system. The model has three channel sections. There is an inflow forcing at the upstream boundary and a water level at the downstream boundary. The decision variables are the flow rates (and by extension, the weir settings) between the channels.

In OpenModelica Connection Editor, the model looks like this:

../../_images/channel_wave_damping_omedit.png

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

 1model Example
 2  // Structures
 3  Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Discharge Discharge;
 4  Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Level Level;
 5  Deltares.ChannelFlow.Hydraulic.Branches.HomotopicTrapezoidal upstream(
 6    theta = theta,
 7    semi_implicit_step_size = step_size,
 8    H_b_up = 15,
 9    H_b_down = 15,
10    bottom_width_up = 50,
11    bottom_width_down = 50,
12    length = 20000,
13    uniform_nominal_depth = 5,
14    friction_coefficient = 35,
15    n_level_nodes = 4,
16    Q_nominal = 100.0
17  );
18  Deltares.ChannelFlow.Hydraulic.Branches.HomotopicTrapezoidal middle(
19    theta = theta,
20    semi_implicit_step_size = step_size,
21    H_b_up = 10,
22    H_b_down = 10,
23    bottom_width_up = 50,
24    bottom_width_down = 50,
25    length = 20000,
26    uniform_nominal_depth = 5,
27    friction_coefficient = 35,
28    n_level_nodes = 4,
29    Q_nominal = 100.0
30  );
31  Deltares.ChannelFlow.Hydraulic.Branches.HomotopicTrapezoidal downstream(
32    theta = theta,
33    semi_implicit_step_size = step_size,
34    H_b_up = 5,
35    H_b_down = 5,
36    bottom_width_up = 50,
37    bottom_width_down = 50,
38    length = 20000,
39    uniform_nominal_depth = 5,
40    friction_coefficient = 35,
41    n_level_nodes = 4,
42    Q_nominal = 100.0
43  );
44  Deltares.ChannelFlow.Hydraulic.Structures.DischargeControlledStructure dam_middle;
45  Deltares.ChannelFlow.Hydraulic.Structures.DischargeControlledStructure dam_upstream;
46
47  // Inputs
48  input Modelica.SIunits.Position Level_H(fixed = true) = Level.H;
49  input Modelica.SIunits.VolumeFlowRate Inflow_Q(fixed = true) = Discharge.Q;
50
51  // Outputs
52  output Modelica.SIunits.Position H_middle = middle.H[middle.n_level_nodes];
53  output Modelica.SIunits.Position H_upstream = upstream.H[upstream.n_level_nodes];
54  output Modelica.SIunits.VolumeFlowRate Q_in = Inflow_Q;
55
56  // Parameters
57  parameter Modelica.SIunits.Duration step_size;
58  parameter Real theta;
59equation
60  connect(dam_middle.HQDown, downstream.HQUp);
61  connect(dam_middle.HQUp, middle.HQDown);
62  connect(dam_upstream.HQDown, middle.HQUp);
63  connect(dam_upstream.HQUp, upstream.HQDown);
64  connect(Discharge.HQ, upstream.HQUp);
65  connect(Level.HQ, downstream.HQDown);
66initial equation
67  downstream.Q[2:downstream.n_level_nodes + 1] = Inflow_Q;
68  middle.Q[2:middle.n_level_nodes + 1] = Inflow_Q;
69  upstream.Q[2:upstream.n_level_nodes + 1] = Inflow_Q;
70end Example;

Note

In order to simulate and show how the PID controllers activate only once the incoming wave has propagated downstream, we will discretize the model in time with a resolution of 5 minutes. With our spatial resolution of 4 level nodes per 20 km reach, this produces a CFL number of approximately 0.4.

For optimization-based control, such a fine temporal resolution is not needed, as the system is able to look ahead and plan corrective measures ahead of time. In this case, CFL numbers of up to 1 or even higher are typically used.

Nevertheless, in order to present a consistent comparison, a 5 minute time step is also used for the optimization example. It is easy to explore the effect of the time step size on the optimization results by changing the value of the step_size class variable.

To run the model with the local control scheme, we make a second model that constrains the flow rates to the weir settings as determined by local PID controller elements:

 1model ExampleLocalControl
 2  extends Example;
 3  // Add PID Controllers that apply local control
 4  PIDController upstream_pid(
 5    state = dam_upstream.HQUp.H,
 6    target_value = 20.0,
 7    P = -200.0,
 8    I = -0.01,
 9    D = 0.0,
10    feed_forward = 100.0,
11    control_action = dam_upstream.Q
12  );
13  PIDController middle_pid(
14    state = dam_middle.HQUp.H,
15    target_value = 15.0,
16    P = -200.0,
17    I = -0.01,
18    D = 0.0,
19    feed_forward = 100.0,
20    control_action = dam_middle.Q
21  );
22  output Modelica.SIunits.VolumeFlowRate Q_dam_upstream = dam_upstream.Q;
23  output Modelica.SIunits.VolumeFlowRate Q_dam_middle = dam_middle.Q;
24end ExampleLocalControl;

The local control model makes use of a PID controller class:

 1model PIDController
 2  input Real state;
 3  parameter Real target_value;
 4  parameter Real P = 1.0;
 5  parameter Real I = 0.0;
 6  parameter Real D = 0.0;
 7  parameter Real feed_forward = 0.0;
 8  output Real control_action;
 9  Real _error;
10  Real error_integral(nominal = 3600);
11equation
12  _error = target_value - state;
13  der(error_integral) = _error;
14  control_action = feed_forward + P * _error + I * error_integral + D * der(_error);
15initial equation
16  error_integral = 0.0;
17end PIDController;

Important

Modellers should take care to set proper values for the initial derivatives, in order to avoid spurious waves at the start of the optimization run. In this example we assume a steady state initial condition, as indicated and enforced by the SteadyStateInitializationMixin in the Python code.

The Optimization Problem

Goals

In this model, we define a TargetLevelGoal to find a requested target level:

14class TargetLevelGoal(Goal):
15    """Really Simple Target Level Goal"""
16
17    def __init__(self, state, target_level):
18        self.function_range = target_level - 5.0, target_level + 5.0
19        self.function_nominal = target_level
20        self.target_min = target_level
21        self.target_max = target_level
22        self.state = state
23
24    def function(self, optimization_problem, ensemble_member):
25        return optimization_problem.state(self.state)
26
27    priority = 1

We will later apply this goal to the upstream and middle channels.

You can read more about the components of goals in the documentation: Multi-objective optimization.

Optimization Problem

We construct the class by declaring it and inheriting the desired parent classes.

30class ExampleOptimization(
31    StepSizeParameterMixin,
32    SteadyStateInitializationMixin,
33    HomotopyMixin,
34    GoalProgrammingMixin,
35    CSVMixin,
36    ModelicaMixin,
37    CollocatedIntegratedOptimizationProblem,
38):

The StepSizeParameterMixin defines the step size parameter and sets the optimization time steps, while the SteadyStateInitializationMixin constrains the initial conditions to be steady-state.

Next, we instantiate the goals. There are two water level goals, applied at the upper and middle channels. The goals are very simple—they just target a specific water level.

41    def path_goals(self):
42        # Add water level goals
43        return [
44            TargetLevelGoal("dam_upstream.HQUp.H", 20.0),
45            TargetLevelGoal("dam_middle.HQUp.H", 15.0),
46        ]

We want to apply these goals to every timestep, so we use the path_goals() method. This is a method that returns a list of the path goals we defined above. Note that with path goals, each timestep is implemented as an independent goal— if we cannot satisfy our min/max on time step A, it will not affect our desire to satisfy the goal at time step B.

For comparison, we also define an optimization problem that uses a local control scheme. This example does not use any goals, as the flow rate is regulated by the PID Controller.

12class ExampleLocalControl(
13    StepSizeParameterMixin,
14    SteadyStateInitializationMixin,
15    HomotopyMixin,
16    CSVMixin,
17    ModelicaMixin,
18    CollocatedIntegratedOptimizationProblem,
19):
20    """Local Control Approach"""
21
22    timeseries_export_basename = "timeseries_export_local_control"

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 classes we just created. We do this for both the local control model and the goal programming model.

50run_optimization_problem(ExampleOptimization)
51run_optimization_problem(ExampleLocalControl)

The Whole Script

All together, all the scripts are as as follows:

step_size_parameter_mixin.py:

 1import numpy as np
 2from rtctools.optimization.optimization_problem import OptimizationProblem
 3
 4
 5class StepSizeParameterMixin(OptimizationProblem):
 6    step_size = 5 * 60  # 5 minutes
 7
 8    def times(self, variable=None):
 9        times = super().times(variable)
10        return np.arange(times[0], times[-1], self.step_size)
11
12    def parameters(self, ensemble_member):
13        p = super().parameters(ensemble_member)
14        p["step_size"] = self.step_size
15        return p

steady_state_initialization_mixin.py:

 1from rtctools.optimization.optimization_problem import OptimizationProblem
 2
 3
 4class SteadyStateInitializationMixin(OptimizationProblem):
 5    def constraints(self, ensemble_member):
 6        c = super().constraints(ensemble_member)
 7        times = self.times()
 8        parameters = self.parameters(ensemble_member)
 9        # Force steady-state initialization at t0 and at t1.
10        for reach in ["upstream", "middle", "downstream"]:
11            for i in range(int(parameters["{}.n_level_nodes".format(reach)]) + 1):
12                state = "{}.Q[{}]".format(reach, i + 1)
13                c.append((self.der_at(state, times[0]), 0, 0))
14        return c

example_local_control.py:

 1from rtctools.optimization.collocated_integrated_optimization_problem import (
 2    CollocatedIntegratedOptimizationProblem,
 3)
 4from rtctools.optimization.csv_mixin import CSVMixin
 5from rtctools.optimization.homotopy_mixin import HomotopyMixin
 6from rtctools.optimization.modelica_mixin import ModelicaMixin
 7from rtctools.util import run_optimization_problem
 8from steady_state_initialization_mixin import SteadyStateInitializationMixin
 9from step_size_parameter_mixin import StepSizeParameterMixin
10
11
12class ExampleLocalControl(
13    StepSizeParameterMixin,
14    SteadyStateInitializationMixin,
15    HomotopyMixin,
16    CSVMixin,
17    ModelicaMixin,
18    CollocatedIntegratedOptimizationProblem,
19):
20    """Local Control Approach"""
21
22    timeseries_export_basename = "timeseries_export_local_control"
23
24
25if __name__ == "__main__":
26    run_optimization_problem(ExampleLocalControl)

example_optimization.py:

 1from example_local_control import ExampleLocalControl
 2from rtctools.optimization.collocated_integrated_optimization_problem import (
 3    CollocatedIntegratedOptimizationProblem,
 4)
 5from rtctools.optimization.csv_mixin import CSVMixin
 6from rtctools.optimization.goal_programming_mixin import Goal, GoalProgrammingMixin
 7from rtctools.optimization.homotopy_mixin import HomotopyMixin
 8from rtctools.optimization.modelica_mixin import ModelicaMixin
 9from rtctools.util import run_optimization_problem
10from steady_state_initialization_mixin import SteadyStateInitializationMixin
11from step_size_parameter_mixin import StepSizeParameterMixin
12
13
14class TargetLevelGoal(Goal):
15    """Really Simple Target Level Goal"""
16
17    def __init__(self, state, target_level):
18        self.function_range = target_level - 5.0, target_level + 5.0
19        self.function_nominal = target_level
20        self.target_min = target_level
21        self.target_max = target_level
22        self.state = state
23
24    def function(self, optimization_problem, ensemble_member):
25        return optimization_problem.state(self.state)
26
27    priority = 1
28
29
30class ExampleOptimization(
31    StepSizeParameterMixin,
32    SteadyStateInitializationMixin,
33    HomotopyMixin,
34    GoalProgrammingMixin,
35    CSVMixin,
36    ModelicaMixin,
37    CollocatedIntegratedOptimizationProblem,
38):
39    """Goal Programming Approach"""
40
41    def path_goals(self):
42        # Add water level goals
43        return [
44            TargetLevelGoal("dam_upstream.HQUp.H", 20.0),
45            TargetLevelGoal("dam_middle.HQUp.H", 15.0),
46        ]
47
48
49# Run
50run_optimization_problem(ExampleOptimization)
51run_optimization_problem(ExampleLocalControl)

Extracting Results

The results from the run are found in output/timeseries_export.csv and output/timeseries_export_local_control.csv. Here are the results when plotted using the python library matplotlib:

(Source code, svg, png)

../../_images/channel_wave_damping.svg

In this example, the PID controller is tuned poorly and ends up amplifying the incoming wave as it propagates downstream. The optimizing controller, in contrast, does not amplify the wave and maintains the target water level throughout the wave event.