...
 
Commits (4)
......@@ -12,7 +12,7 @@ model Example
length = 20000,
uniform_nominal_depth = 5,
friction_coefficient = 35,
n_level_nodes = 4,
n_level_nodes = 21,
Q(each nominal = 100.0)
) annotation(Placement(visible = true, transformation(origin = {-60, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
Deltares.ChannelFlow.Hydraulic.Branches.HomotopicTrapezoidal middle(
......@@ -25,7 +25,7 @@ model Example
length = 20000,
uniform_nominal_depth = 5,
friction_coefficient = 35,
n_level_nodes = 4,
n_level_nodes = 21,
Q(each nominal = 100.0)
) annotation(Placement(visible = true, transformation(origin = {0, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
Deltares.ChannelFlow.Hydraulic.Branches.HomotopicTrapezoidal downstream(
......@@ -38,7 +38,7 @@ model Example
length = 20000,
uniform_nominal_depth = 5,
friction_coefficient = 35,
n_level_nodes = 4,
n_level_nodes = 21,
Q(each nominal = 100.0)
) annotation(Placement(visible = true, transformation(origin = {58, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
Deltares.ChannelFlow.Hydraulic.Structures.DischargeControlledStructure dam_middle annotation(Placement(visible = true, transformation(origin = {30, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
......@@ -64,7 +64,12 @@ equation
connect(Discharge.HQ, upstream.HQUp) annotation(Line(points = {{-82, 0}, {-68, 0}}, color = {0, 0, 255}));
connect(Level.HQ, downstream.HQDown) annotation(Line(points = {{66, 0}, {82, 0}, {82, 0}, {82, 0}}, color = {0, 0, 255}));
initial equation
downstream.Q[2:downstream.n_level_nodes + 1] = Inflow_Q;
middle.Q[2:middle.n_level_nodes + 1] = Inflow_Q;
upstream.Q[2:upstream.n_level_nodes + 1] = Inflow_Q;
//der(downstream.H[1:upstream.n_level_nodes - 1]) = 0;//fill(0.0, upstream.n_level_nodes - 1);
//der(middle.H[1:middle.n_level_nodes - 1]) = 0; //fill(0.0, middle.n_level_nodes - 1);
//der(upstream.H[1:downstream.n_level_nodes - 1]) = 0;// fill(0.0, downstream.n_level_nodes - 1);
//downstream.Q[2:downstream.n_level_nodes + 1] = Inflow_Q;
//middle.Q[2:middle.n_level_nodes + 1] = Inflow_Q;
//upstream.Q[2:upstream.n_level_nodes + 1] = Inflow_Q;
// Start on set point
//downstream.H[downstream.n_level_nodes + 0] = 10;
end Example;
model ExampleHybrid
extends Example;
PIDController middle_pid(
state = dam_middle.HQUp.H,
target_value = 15.0,
P = -200.0,
I = -0.01,
D = 0.0,
feed_forward = 100.0,
control_action = dam_middle.Q
);
output Modelica.SIunits.VolumeFlowRate Q_dam_middle = dam_middle.Q;
input Modelica.SIunits.VolumeFlowRate Q_dam_upstream(fixed = false, min = 0.0, max = 1000.0, nominal = 100.0) = dam_upstream.Q;
initial equation
//upstream.H[4] = 20.0;
//middle.H[middle.n_level_nodes] = 15;
upstream.H[upstream.n_level_nodes] = 20;
end ExampleHybrid;
model ExampleOptimization
extends Example;
input Modelica.SIunits.VolumeFlowRate Q_dam_middle(fixed = false, min = 0.0, max = 1000.0, nominal = 100.0) = dam_middle.Q;
input Modelica.SIunits.VolumeFlowRate Q_dam_upstream(fixed = false, min = 0.0, max = 1000.0, nominal = 100.0) = dam_upstream.Q;
input Modelica.SIunits.VolumeFlowRate Q_dam_middle(fixed = false, min = 0.0, max = 1000.0, nominal = 100.0) = dam_middle.Q;
initial equation
middle.H[middle.n_level_nodes] = 15;
upstream.H[upstream.n_level_nodes] = 20;
end Example;
......@@ -7,7 +7,7 @@ model PIDController
parameter Real feed_forward = 0.0;
output Real control_action;
Real _error;
Real _error_integral(nominal = 3600);
Real _error_integral(nominal = 1e6);
equation
_error = target_value - state;
der(_error_integral) = _error;
......
import numpy as np
import casadi as ca
from example_local_control import ExampleLocalControl
from rtctools.optimization.collocated_integrated_optimization_problem import (
CollocatedIntegratedOptimizationProblem,
)
from rtctools.optimization.csv_mixin import CSVMixin
from rtctools.optimization.goal_programming_mixin import Goal, GoalProgrammingMixin
from rtctools.optimization.homotopy_mixin import HomotopyMixin
from rtctools.optimization.modelica_mixin import ModelicaMixin
from rtctools.util import run_optimization_problem
from steady_state_initialization_mixin import SteadyStateInitializationMixin
from step_size_parameter_mixin import StepSizeParameterMixin
class TargetLevelGoal(Goal):
"""Really Simple Target Level Goal"""
def __init__(self, state, target_level, margin):
self.function_range = target_level - 10, target_level + 10
self.function_nominal = target_level
self.target_min = target_level - margin / 2
self.target_max = target_level + margin / 2
self.state = state
def function(self, optimization_problem, ensemble_member):
return optimization_problem.state(self.state)
priority = 1
class ExampleHybrid(
StepSizeParameterMixin,
SteadyStateInitializationMixin,
HomotopyMixin,
GoalProgrammingMixin,
CSVMixin,
ModelicaMixin,
CollocatedIntegratedOptimizationProblem,
):
"""Goal Programming Approach"""
def path_goals(self):
# Add water level goals
goals = [
TargetLevelGoal("dam_upstream.HQUp.H", 20.0, 0),
]
return goals
def __solver_options(self):
options = super().solver_options()
options['ipopt.linear_solver'] = 'ma86'
options['ipopt.nlp_scaling_method'] = 'none'
options['ipopt.jac_c_constant'] = 'no'
options['ipopt.jac_d_constant'] = 'yes'
options['expand'] = True #TODO
return options
def __goal_programming_options(self):
options = super().goal_programming_options()
options['scale_by_problem_size'] = False
options['keep_soft_constraints'] = False
return options
# Run
run_optimization_problem(ExampleHybrid)
#run_optimization_problem(ExampleLocalControl)
import numpy as np
import casadi as ca
from example_local_control import ExampleLocalControl
from rtctools.optimization.collocated_integrated_optimization_problem import (
CollocatedIntegratedOptimizationProblem,
)
from rtctools.optimization.csv_mixin import CSVMixin
from rtctools.optimization.goal_programming_mixin import Goal, GoalProgrammingMixin
from rtctools.optimization.homotopy_mixin import HomotopyMixin
from rtctools.optimization.modelica_mixin import ModelicaMixin
from rtctools.util import run_optimization_problem
from steady_state_initialization_mixin import SteadyStateInitializationMixin
from step_size_parameter_mixin import StepSizeParameterMixin
class TargetLevelGoal(Goal):
"""Really Simple Target Level Goal"""
def __init__(self, state, target_level, margin, priority):
self.function_range = target_level - 10, target_level + 10
self.function_nominal = 10
self.target_min = target_level - margin / 2
self.target_max = target_level + margin / 2
self.state = state
self.priority = priority
def function(self, optimization_problem, ensemble_member):
return optimization_problem.state(self.state)
order = 2
class TargetMaximumDischargeGoal(Goal):
"""Really Simple Target Level Goal"""
def __init__(self, state, target_level, margin, priority):
self.function_range = 0, 1000
self.function_nominal = 100
self.target_max = target_level + margin / 2
self.state = state
self.priority = priority
def function(self, optimization_problem, ensemble_member):
return optimization_problem.state(self.state)
order = 2
class SmoothingGoal(Goal):
"""Smoothing Goal"""
def __init__(self, state, nominal, priority):
self.function_nominal = nominal
self.state = state
self.priority = priority
def function(self, optimization_problem, ensemble_member):
return optimization_problem.der(self.state)
order = 2
class MinAmplitudeGoal(Goal):
"""Amplitude Minimization Goal"""
def __init__(self, variable, priority):
self.function_nominal = 100.0
self.variable = variable
self.priority = priority
def function(self, optimization_problem, ensemble_member):
return optimization_problem.extra_variable(self.variable)
order = 2
class ExampleOptimization(
StepSizeParameterMixin,
SteadyStateInitializationMixin,
HomotopyMixin,
GoalProgrammingMixin,
CSVMixin,
ModelicaMixin,
CollocatedIntegratedOptimizationProblem,
):
"""Goal Programming Approach"""
def pre(self):
super().pre()
self._a_u = ca.MX.sym('a_u')
self._a_m = ca.MX.sym('a_m')
@property
def extra_variables(self):
return super().extra_variables + [self._a_u, self._a_m]
def bounds(self):
bounds = super().bounds()
bounds['a_u'] = (0.0, 1000.0)
bounds['a_m'] = (0.0, 1000.0)
return bounds
def variable_nominal(self, variable):
if variable in set(['a_u', 'a_m']):
return 100.0
else:
return super().variable_nominal(variable)
def path_constraints(self, ensemble_member):
path_constraints = super().path_constraints(ensemble_member)
path_constraints.append((self.state("dam_upstream.HQUp.Q") - self._a_u, -np.inf, 0.0))
path_constraints.append((self.state("dam_middle.HQUp.Q") - self._a_m, -np.inf, 0.0))
return path_constraints
def goals(self):
goals = [
MinAmplitudeGoal('a_u', 2),
MinAmplitudeGoal('a_m', 2),
]
return goals
def path_goals(self):
path_goals = [
TargetLevelGoal("dam_upstream.HQUp.H", 20.0, 1.0, 1),
TargetLevelGoal("dam_middle.HQUp.H", 15.0, 1.0, 1),
TargetLevelGoal("dam_upstream.HQUp.H", 20.0, 0.0, 2),
TargetLevelGoal("dam_middle.HQUp.H", 15.0, 0.0, 2),
#TargetMaximumDischargeGoal("dam_upstream.HQUp.Q", 100, 0.0, 2), # instead of MinAmplitudeGoals above
#TargetMaximumDischargeGoal("dam_middle.HQUp.Q", 100, 0.0, 2),
SmoothingGoal("dam_upstream.HQUp.Q", 1e-2, 4),
SmoothingGoal("dam_middle.HQUp.Q", 1e-2, 4),
]
# TODO idea: spmile cutoff goal
return path_goals
def solver_options(self):
options = super().solver_options()
options['ipopt.linear_solver'] = 'ma86'
options['ipopt.nlp_scaling_method'] = 'none'
options['ipopt.jac_c_constant'] = 'no'
options['ipopt.jac_d_constant'] = 'yes'
options['expand'] = True # TODO port to other examples.
return options
def goal_programming_options(self):
options = super().goal_programming_options()
options['keep_soft_constraints'] = True
options['scale_by_problem_size'] = False
return options
# Run
run_optimization_problem(ExampleOptimization)
#run_optimization_problem(ExampleLocalControl)
import numpy as np
import casadi as ca
from example_local_control import ExampleLocalControl
from rtctools.optimization.collocated_integrated_optimization_problem import (
......@@ -17,17 +20,62 @@ from step_size_parameter_mixin import StepSizeParameterMixin
class TargetLevelGoal(Goal):
"""Really Simple Target Level Goal"""
def __init__(self, state, target_level):
self.function_range = target_level - 5.0, target_level + 5.0
self.function_nominal = target_level
self.target_min = target_level
self.target_max = target_level
def __init__(self, state, target_level, margin, priority):
self.function_range = target_level - 10, target_level + 10
self.function_nominal = 10
self.target_min = target_level - margin / 2
self.target_max = target_level + margin / 2
self.state = state
self.priority = priority
def function(self, optimization_problem, ensemble_member):
return optimization_problem.state(self.state)
order = 2
class TargetMaximumDischargeGoal(Goal):
"""Really Simple Target Level Goal"""
def __init__(self, state, target_level, margin, priority):
self.function_range = 0, 1000
self.function_nominal = 100
self.target_max = target_level + margin / 2
self.state = state
self.priority = priority
def function(self, optimization_problem, ensemble_member):
return optimization_problem.state(self.state)
priority = 1
order = 2
class SmoothingGoal(Goal):
"""Smoothing Goal"""
def __init__(self, state, nominal, priority):
self.function_nominal = nominal
self.state = state
self.priority = priority
def function(self, optimization_problem, ensemble_member):
return optimization_problem.der(self.state)
order = 2
class MinAmplitudeGoal(Goal):
"""Amplitude Minimization Goal"""
def __init__(self, variable, priority):
self.function_nominal = 100.0
self.variable = variable
self.priority = priority
def function(self, optimization_problem, ensemble_member):
return optimization_problem.extra_variable(self.variable)
order = 2
class ExampleOptimization(
......@@ -41,14 +89,82 @@ class ExampleOptimization(
):
"""Goal Programming Approach"""
def pre(self):
super().pre()
self._a_u = ca.MX.sym('a_u')
self._a_m = ca.MX.sym('a_m')
@property
def extra_variables(self):
return super().extra_variables #+ [self._a_u, self._a_m]
def bounds(self):
bounds = super().bounds()
bounds['a_u'] = (0.0, 1000.0)
bounds['a_m'] = (0.0, 1000.0)
return bounds
def constant_inputs(self, ensemble_member):
constant_inputs = super().constant_inputs(ensemble_member)
#constant_inputs['Inflow_Q'].values.fill(500.0)
#constant_inputs['Inflow_Q'].values[0] = 100.0
#constant_inputs['Inflow_Q'].values.fill(100.0)
#constant_inputs['Inflow_Q'].values[0] = 500.0
return constant_inputs
def variable_nominal(self, variable):
if variable in set(['a_u', 'a_m']):
return 100.0
else:
return super().variable_nominal(variable)
def path_constraints(self, ensemble_member):
path_constraints = super().path_constraints(ensemble_member)
#path_constraints.append((self.state("dam_upstream.HQUp.Q") - self._a_u, -np.inf, 0.0))
#path_constraints.append((self.state("dam_middle.HQUp.Q") - self._a_m, -np.inf, 0.0))
return path_constraints
def goals(self):
goals = [
#MinAmplitudeGoal('a_u', 2),
#MinAmplitudeGoal('a_m', 2),
]
return goals
def path_goals(self):
# Add water level goals
return [
TargetLevelGoal("dam_upstream.HQUp.H", 20.0),
TargetLevelGoal("dam_middle.HQUp.H", 15.0),
path_goals = [
#TargetLevelGoal("dam_upstream.HQUp.H", 20.0, 1.0, 1),
#TargetLevelGoal("dam_middle.HQUp.H", 15.0, 1.0, 1),
TargetLevelGoal("dam_upstream.HQUp.H", 20.0, 0.0, 2),
TargetLevelGoal("dam_middle.HQUp.H", 15.0, 0.0, 2),
#TargetMaximumDischargeGoal("dam_upstream.HQUp.Q", 100, 0.0, 2),
#TargetMaximumDischargeGoal("dam_middle.HQUp.Q", 100, 0.0, 2),
#SmoothingGoal("dam_upstream.HQUp.Q", 1e-2, 4),
#SmoothingGoal("dam_middle.HQUp.Q", 1e-2, 4),
]
# TODO idea: spmile cutoff goal
return path_goals
def solver_options(self):
options = super().solver_options()
options['ipopt.linear_solver'] = 'ma86'
options['ipopt.nlp_scaling_method'] = 'none'
options['ipopt.jac_c_constant'] = 'no'
options['ipopt.jac_d_constant'] = 'yes'
options['expand'] = True # TODO port to other examples.
return options
def goal_programming_options(self):
options = super().goal_programming_options()
options['keep_soft_constraints'] = True
options['scale_by_problem_size'] = False
return options
# Run
run_optimization_problem(ExampleOptimization)
run_optimization_problem(ExampleLocalControl)
#run_optimization_problem(ExampleLocalControl)
......@@ -7,10 +7,17 @@ class SteadyStateInitializationMixin(OptimizationProblem):
times = self.times()
parameters = self.parameters(ensemble_member)
# Force steady-state initialization at t0 and at t1.
# TODO move to initial equation section.
n_level_nodes = 21
for reach in ['upstream', 'middle', 'downstream']:
for i in range(int(parameters[f'{reach}.n_level_nodes']) + 1):
for i in range(n_level_nodes + 1):
state = f'{reach}.Q[{i + 1}]'
c.append(
(self.der_at(state, times[0]), 0, 0)
)
for i in range(n_level_nodes):
state = f'{reach}.H[{i + 1}]'
c.append(
(self.der_at(state, times[0]), 0, 0)
)
return c
......@@ -4,7 +4,7 @@ from rtctools.optimization.optimization_problem import OptimizationProblem
class StepSizeParameterMixin(OptimizationProblem):
step_size = 5 * 60 # 5 minutes
step_size = 15 * 60 # 15 minutes
def times(self, variable=None):
times = super().times(variable)
......@@ -14,3 +14,14 @@ class StepSizeParameterMixin(OptimizationProblem):
p = super().parameters(ensemble_member)
p['step_size'] = self.step_size
return p
def solver_options(self):
o = super().solver_options()
o['expand'] = True
o['ipopt.linear_solver'] = 'ma86'
return o
def compiler_options(self):
o = super().compiler_options()
o['replace_parameter_values'] = True
return o
......@@ -512,6 +512,9 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass=ABC
constant_parameters,
constant_parameter_values)
# Collect extra variable symbols
symbolic_extra_variables = ca.vertcat(*self.extra_variables)
# Aggregate ensemble data
ensemble_aggregate = {}
ensemble_aggregate["parameters"] = ca.horzcat(*[nullvertcat(*l) for l in ensemble_parameter_values])
......@@ -713,7 +716,8 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass=ABC
*integrated_variables, *collocated_variables, *integrated_derivatives,
*collocated_derivatives, *self.dae_variables['constant_inputs'],
*self.dae_variables['time'], *self.path_variables,
*self.__extra_constant_inputs)]
*self.__extra_constant_inputs),
symbolic_extra_variables]
# Initialize a Function for the path objective
# Note that we assume that the path objective expression is the same for all ensemble members
......@@ -769,6 +773,7 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass=ABC
collocation_time_0 = accumulated_U[offset + 0]
collocation_time_1 = accumulated_U[offset + 1]
path_variables_1 = accumulated_U[offset + 2:offset + 2 + len(self.path_variables)]
extra_variables = ca.vertcat(*[self.extra_variable(var.name()) for var in self.extra_variables])
extra_constant_inputs_1 = accumulated_U[offset + 2 + len(self.path_variables):]
# Approximate derivatives using backwards finite differences
......@@ -856,7 +861,8 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass=ABC
constant_inputs_1,
collocation_time_1 - t0,
path_variables_1,
extra_constant_inputs_1)]
extra_constant_inputs_1),
symbolic_extra_variables]
accumulated_Y.extend(path_objective_function.call(
self.__func_inputs_implicit, False, True))
......@@ -868,7 +874,9 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass=ABC
self.__func_inputs_implicit, False, True))
# Save the accumulated inputs such that can be used later in map_path_expression()
self.__func_accumulated_inputs = (accumulated_X, accumulated_U, symbolic_parameters)
self.__func_accumulated_inputs = (
accumulated_X, accumulated_U,
ca.veccat(symbolic_parameters, symbolic_extra_variables))
# Use map/mapaccum to capture integration and collocation constraint generation over the
# entire time horizon with one symbolic operation. This saves a lot of memory.
......@@ -1090,11 +1098,13 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass=ABC
# Save these inputs such that can be used later in map_path_expression()
self.__func_mapped_inputs.append(
(accumulation_X0, accumulation_U, ca.repmat(parameters, 1, n_collocation_times - 1)))
(accumulation_X0, accumulation_U,
ca.repmat(ca.vertcat(parameters, extra_variables), 1, n_collocation_times - 1)))
self.__func_initial_inputs.append([parameters, ca.vertcat(
initial_state, initial_derivatives, initial_constant_inputs, 0.0,
initial_path_variables, initial_extra_constant_inputs)])
initial_path_variables, initial_extra_constant_inputs),
extra_variables])
integrators_and_collocation_and_path_constraints = accumulation(
*self.__func_mapped_inputs[ensemble_member])
......@@ -1200,7 +1210,7 @@ class CollocatedIntegratedOptimizationProblem(OptimizationProblem, metaclass=ABC
initial_constant_inputs,
0.0,
[self.variable_nominal(var.name()) for var in self.path_variables],
initial_extra_constant_inputs)])
initial_extra_constant_inputs), extra_variables])
for i in range(len(delayed_feedback_expressions)):
expression = delayed_feedback_expressions[i]
......