-
Notifications
You must be signed in to change notification settings - Fork 2
Home
PARODIS is a MATLAB framework for Pareto Optimal (scenario- based economic) Model Predictive Control for Distributed Systems. The main strength of PARODIS is, that the underlying optimization problem can be formulated completely symbolically and parametrically. PARODIS uses the YALMIP library by Johan Löfberg for the symbolic representation of optimization problems and for interfacing the numerical solvers.
PARODIS efficiently manages the symbolic representation and uses the optimizer
feature of YALMIP to pre-compile the symbolically defined optimization problem, which makes PARODIS rather fast.
In PARODIS, a simulation consists of multiple agents, where each agent represents it's own (sub)system with it's own model and corresponding controller. At the beginning of each global simulation step, these agents can interact and exchange information, e.g. for consensus algorithms or hierarchical information exchange. In PARODIS, this interaction is called a negotiation. After this negotiation, the actual time steps of the agent are executed.
The model representation, which is the basis of each agent's system, is a symbolic state space model, with a state
In PARODIS, the controllers are by default (economic) MPC controllers. A control strategy is defined by a set of constraints, symbolic parameters and a cost function that is to be minimized. The controller manages these constraints, parameters and cost functions. Parameters are expressed symbolically and defined a data source, which is used to replace parameters dynamically at runtime.
Cost functions are classes, that can yield scalar symbolic expression for costs over the horizon based on a symbolic model and parameter definition. Cost functions can introduce additional constraints and slack variables for an enhanced (economic) optimization problem formulation.
Since PARODIS offers the functionality of Pareto optimal MPC, multiple cost functions can be added to the controller. These cost functions can be weighted against each other; either manually, as is the classical approach in (economic) MPC, or dynamically by Pareto evaluation.
PARODIS also allows the definition of evaluation functions, which are called to evaluate both the current prediction horizon as well as a performed simulation step. This allows for a live evaluation of the simulation, making a posteriori analysis easier and faster.
Furthermore, PARODIS has a powerful and easy-to-use plotting interface, with which all important information, including the aforementioned evaluation functions, can be plotted - both live, i.e. while the simulation is running, and after the simulation is finished.
- Getting Started
PARODIS is compatible with MATLAB 2020a down to MATLAB 2018a. 2018a is the lowest version PARODIS has been tested with.
The framework does not require any special toolboxes.
PARODIS does strongly depend on YALMIP by Johan Löfberg and was developed with YALMIP version 20200116. In order for PARODIS to work, YALMIP must be installed and avaiable in your MATLAB path. See the YALMIP project website for install instructions and tips and tricks on how to efficiently formulate optimization problems.
In order to solve optimisation problems, YALMIP interfaces a wide variety of solvers. Please note, that no solvers are provided with YALMIP or PARODIS, and they must be installed an configured seperately.
Installing and setting up PARODIS is quite straight forward. First, install YALMIP according to its install instructions and install any solvers you may need.
To install PARODIS, download the latest release or clone this repository to your desired destination. Then, open MATLAB and run the install script. This script will add the correct directories to your MATLAB path and check if YALMIP is installed.
cd YourPARODISdirectoryGoesHere
install
In the language of PARODIS, an Ansatz is the combination of a model with an (MPC) controller for a certain control problem, given a certain prediction horizon and a number of scenarios that should be factored in.
A model is a symbolic representation of the state trajectory over the horizon, together with a symbolic representation of the input u and the disturbances d. A model is defined by a set of ordinary difference equations for the states that yield
If the argument implicitPrediction
is set to false, the state trajectory
To create a model with this definition goes as follows:
T_s = [5 5 10 10 60]; % horizon with non equidistant timesteps
numScenarios = 10; % number of scenarios to be evaluated
implicitPrediction = true;
model = createModel( @my_model_fun, T_s, numScenarios, implicitPrediction );
where [ode, n_x, n_u, n_d] = my_model_fun( T_s_n )
returns an ODE and the number of states/inputs/disturbances respectively for a given sample time T_s_n
.
function [ode, n_x, n_u, n_d] = my_model_fun(T_s)
ode = @(x, u, d)( [ x(1)^2 + u(1) + sin(d(1)); ...
x(2) + u(2) + cos(d(1)) ] );
n_x = 2;
n_u = 2;
n_d = 2;
end
Note: For linear systems of form
function [ode, n_x, n_u, n_d, linearRepresentation] = test_model(T_s)
A = diag([ 0.1 0.05 1.01 ]);
B = [1 0; 0 1; 1 1];
S = diag([2 2 2]);
% always required!
ode = @(x, u, d)( A*x + B*u + S*d );
n_x = 3;
n_u = 2;
n_d = 3;
% additional linear representation
linearRepresentation = struct;
linearRepresentation.A = A;
linearRepresentation.B = B;
linearRepresentation.S = S;
end
See the Models documentation for more details.
Now, with this, an ansatz can be formulated:
function [model, controller, T_s, numScenarios] = create_my_ansatz(implicitPrediction)
T_s = [5 5 10 10 60]; % horizon with non equidistant timesteps
numScenarios = 10; % number of scenarios to be evaluated
model = createModel( @my_model_fun, T_s, numScenarios, implicitPrediction );
controller = SymbolicController( numScenarios );
myParam = controller.addParam('myParam', ... );
controller.addBoxConstraint("x", ... );
controller.addDeltaConstraint("du", ... );
for s=1:numScenarios
controller.addConstraint( (0 <= model.x{s}(1, :) - model.x{s}(2, :) <= myParam ):'keep state 1 and 2 close' );
end
controller.addCostFunction("energy_costs", EnergyCostFunction);
controller.addCostFunction("co2_emissions", CO2EmissionsCostFunction);
controller.predDisturbanceSource = @my_disturbance_source;
controller.realDisturbanceSource = './data/real_disturbances.csv';
Note that it is not necessary to wrap everything into an ansatz function, but it is very much recommended, especially if you have more than one ansatz you want to test out with your agent or if you are dealing with multiple agents.
An important expression in PARODIS language is the scenario. A scenario
This definition allows the implementation of scenario-based approaches, where an optimal input is to be found and applied to the real system, by incorporating the ensemble of scenarios into the decision (i.e. optimization) process. An example could be the search for an input that minimizes the deviation of the the in a room from some setpoint, averaged over all scenarios, where each scenario presents a different possible outside temperature trajectory.
In PARODIS, you can freely choose how to deal with scenarios in cost functions and evaluation functions, as well as how to incorporate scenario-dependent constraints. Keep in mind that, box and delta constraints (see further down) will automatically be applied equally to all scenarios.
PARODIS offers three different controllers
-
SymbolicController
, where the optimization problem is compiled initially as completely parametric YALMIPoptimizer
object, i.e. the problem remains fully parametrized at all times -
ExplicitController
where the optimization problem is recreated at each time step, parameters are always directly replaced by their values and YALMIP'soptimize()
function is used -
ParetoController
Formulates the MPC problem as a Pareto optimization problem. In order to do that, it recreates the basic optimization at each time step like theExplicitController
and then creates a YALMIPoptimizer
to generate Pareto frontiers efficiently
The SymbolicController
is best suited for smaller problems with few scenarios (if any), as the optimizer
parametrization does not scale well. The ExplicitController
on the other hand scales very well, also for large problems, but is comparatively slow for small problems. It is also better suited for debugging, as the created model is more transparent than the compiled optimizer
object.
Initially, you should try using the SymbolicController
if you don't want to use Pareto optimization, and if model compilation/evaluation is slow, switch to the ExplicitController
. Also the ExplicitController
may be better suited while implementing a new problem, as it is more transparent for debugging.
Regarding Pareto optimization and using the ParetoController
, refer to Pareto optimization and the documentation of the controller itself ParetoController
.
PARODIS works with a symbolic representation of the MPC optimization problem. MPC problems often depend on time variant parameters in the cost functions. To support this, you can add parameters to your problem. These parameters can be accessed by their name and used in constraint, cost functions, and evaluation functions.
paramExpression = controller.addParam( name, dimensions, source, scenario_dependent )
where dimensions
is a 2x1 vector [nrows ncols]
passed to sdpvar
. name
must be a valid fieldname for a struct.
The source
can be either a constant matrix, a function handle or a CSV filename. In case of the latter, it will be treated like a disturbance prediction and always retrieve the entire horizon of data from the sourc. See Sources for a detailed description of data sources.
If source
is a function handle, it should return a cell array with numScenarios
entries, where each cell is the parameter value for the given scenario. A parameter source function should have the following signature:
parameterValue = source_function(T, callingAgent, agents, numScenarios);
Where T
is a vector of length length(T_s)
and contains the horizon time steps in simulation time. callingAgent
is the Agent instance of the agent calling the source function, agents
a struct with all named agents of the current simulation.
The argument scenario_dependent
expresses whether the the parameter differs between scenarios or not. This flag should be set properly, as it can increase performance.
There are three options for adding constraints.
- Box constraints are constraint on the states or the inputs that will be set for all scenarios equally and over the entire length of the horizon
controller.addBoxConstraint( variable = "x"/"u", indexes, lb, ub);
- Delta constraints are constraints on the difference
$\Delta x(n|k) = x(n+1|k) - x(n|k)$ (or$u$ respectively), again set equally for all scenarios and the entire horizon. For the delta constraints on the state,n = 0, ..., N_pred
applies. On the input,n = 0, ..., N_pred - 1
respectively. Note, that fordu(0|k) = u(0|k) - u(k-1)
the inputu(k-1)
of the previous time step is automatically provided by PARODIS. At simulation startk = 0
, PARODIS will either assumeu(k-1)
to be zero, or use an initial value provided to theAgent
constructor (see Agent). The argumentdT
defines the time difference to which the delta constraint's lb/ub relate. This is necessary, since the time steps within the horizon are not necessarily equidistant.
controller.addDeltaConstraint( variable = "dx"/"du", indexes, lb, ub, dT)
- Free constraint expressions, that can express any valid constraint
controller.addConstraint( constraint );
**Note!: ** This way of adding a free constraint expression should be avoided when using the ExplicitController, as it requires applying YALMIP's replace
function to the constraint at each time step, to replace parameters with their values, which is very, very slow.
Alternatively, a free constraint expressions can also be defined implicitly, by providing a function handle, which takes the three input arguments model
, parameters
and slacks
and returns a constraint expression when called:
controller.addConstraint( @(model, parameters, slacks, s)( constraint ) )
Example:
controller.addConstraint( @(model, parameters, slacks, s)( 0 <= model.x{s} <= parameters.param{s} ) );
If you don't want to define the constraint for all scenarios s
but explicitly for some scenario (or when you don't use a scenario-based approach), you can omit the fourth argument:
controller.addConstraint( @(model, parameters, slacks)( 0 <= model.x{1} <= parameters.param{1} ) );
If the upper and/or lower bound on addBoxConstraint
or addDeltaConstraint
is a parameter, the constraint will be set equally over the horizon, and coupled by scenarios. Consider this example for clarification:
param = controller.addParam( 'param', [4 1], @my_param_source, true );
controller.addBoxConstraint("x", 1:4, 0, param);
% is equivalent to:
for s=1:length(param)
controller.addConstraint( 0 <= model.x{s} <= param{s} );
end
For PARODIS, cost functions are subclasses of the abstract CostFunction
class. This way, the cost functions can introduce their own slack variables and constraint expressions. See the documentation page Cost Functions for more detailed information.
A simple example cost function class could look like this:
classdef ExampleCostFunction < CostFunction
methods (Static)
function [slacks] = getSlacks(model, agent, params)
slacks = struct;
end
function [constraints] = getConstraints(model, agent, slacks, params)
constraints = [];
end
function expr = buildExpression(x, u, d, params, Ns, slacks)
expr = 0;
for s=1:Ns
expr = expr + ExampleCostFunction.buildExpressionSingleScen(x{s}, u, d{s}, extractScenario(params, s), slacks);
end
end
function [exprSingle] = buildExpressionSingleScen(x, u, d, params, slacks)
exprSingle = sum(x(1, :));
end
function [horizon] = evaluateHorizon(x, u, d, params, slacks)
Ns = length(x);
horizon = x{1}(1, 1:end-1);
for s=2:Ns
horizon = horizon + x{s}(1, 1:end-1);
end
end
end
end
In some cases you may want to define slack variables when defining model constraints, or you may want to introduce slack variables that are used across multiple cost functions. In this case, these variables can be defined by adding a shared slack variable. This can be done by using
slack = controller.addSharedSlack(name, dimensions, [full = false])
where dimensions
is a 2x1 vector [nrows ncols]
passed to sdpvar
. If nrows = ncols
, YALMIP will automatically assume that a symmetric matrix is used. To instead mark the matrix as being non-symmetric, set the optional argument full = true
. name
must be a valid fieldname for a struct.
In PARODIS, two disturbance sources are to be defined:
-
controller.predDisturbanceSource
, from which the disturbance predictions over the horizon are retrieved -
controller.realDisturbanceSource
, from which the actual disturbance that is to applied to the system is retrieved
Both sources can be either a function handle or a path to a CSV file. In case of a CSV file, the disturbances will be interpolated from the time series entries in the CSV file. Note, that the CSV file must contain enough entries to cover the simulation time and the last horizon, since disturbances are interpolated from the data, but not extrapolated! See Sources for a detailed description of data sources.
In case of the source being a function handle, it behaves the same way as a function handle source for a parameter. The function must return a cell array with numScenarios
entries, where each cell is the parameter value for the given scenario. A source function should have the following signature:
data = source_function(T, callingAgent, agents, numScenarios);
Where T
is a vector of length length(T_s)
and contains the horizon time steps in simulation time. callingAgent
is the Agent instance of the agent calling the source function, agents
a struct with all named agents of the current simulation.
There are two ways to create an agent:
- Directly instantiating an Agent object, providing a name, a model struct, a Controller instance and a value for the initial state of the system:
% horizon time steps
T_s = [ .. .. .. ]
agent = Agent(name, model, controller, T_s, x0, uPrev_0 = zeros(n_u, 1));
Optionally, you can provide an initial value for the input u
as well, which will is considered at k = 0
for any constraints on the difference du(0|k) = u(0|k) - u(k-1)
. By default, this initial value is set to the zero vector.
The first time step in the horizon time steps defines the agent's simulation time step.
After creating an agent, you should configure it to your needs, by setting the configuration and adding evaluation functions accordingly, as well as setting the callbacks for state measurement and negotiation.
[model, controller, T_s, numScenarios] = create_my_ansatz();
agent = Agent('my_agent', model, controller, T_s, x0);
% changing configuration parameters
agent.config.solver = 'quadprog';
agent.config.debugMode = true;
% setting callbacks
agent.callbackMeasurement = @my_measurement_callback;
agent.callbackNegotiation = @my_negotiation_callback;
agent.addEvalFunction( 'electricity_costs', @eval_function_for_electricity_costs );
agent.addEvalFunction( 'get_slack', @eval_function_get_slack_variable );
Evaluation functions are named functions, that will be calculated automatically during the simulation before executing a simulation step and afterwards. They can be plotted the same way as cost functions or state trajectories. These functions are given four parameters:
-
agent
the agent calling the function -
simulation
the simulation instance from which it was called -
predict = true/false
true if the function should return a vector with values over the horizon, false if for the realised timestep -
scenario
the scenario which shall be evaluated. IfaddEvalFunction( name, handle, scenarioDependent = false)
is used,scenario = 1
always applies
An evaluation function should return either a scalar value (if predict = false
) or a vector of length length(T_s)
Since these functions are given direct access to the agent and the simulation, they can be used to evaluate any desired data across multiple agents.
Evaluation functions must not necessarily make sense to be defined for both values of predict
. In that case, they should return NaN
or NaN(1, length(T_s))
respectively.
As an example, here an implementation of an evaluation function that returns the value of a slack variable:
function [value] = eval_function_get_slack_variable(agent, simulation, predict, scenario)
if predict
value = abs( agent.status.slackVariables.my_slack );
else
% slack variable is only relevant for horizon
value = NaN;
end
end
Evaluation functions are added directly to the agent
agent.addEvalFunction(name, function_handle[, scenarioDependent = false])
If the evaluation function is scenario dependent, set the corresponding argument to true. Otherwise, PARODIS will only evaluate the eval function for the first scenario (and use this value for all scenarios).
Simulations are configured and represented using an instance of the Simulation class. Running a simulation is pretty straight forward. Instantiating an instance requires two arguments (and one optional argument)
sim = Simulation(name, simulationTime, [addTimestamp = false])
The name of the simulation will define the name of directory in which the simulation results (agent histories, plots, logs, RNG seeds) will be stored.
The simulationTime
gives the time until which the simulation should be run (in the timescale of the agents).
The optional argument addTimestamp
can be set to true
, which adds the date to the simulation name as result directory, or to 'seconds'
, in that case the date and the current time will be added.
After instantiating, the simulation configuration can be adjusted, i.e. enabling/disabling (live) plots, enabled/disabling storing of plots or history.
To run a simulation, the agents that should be simulated have to created and added using addAgent
or addAgents
. Analogously, the figures that should be plotted are to be created and added using addPlot
or addPlots
.
If the agents support negotiation, a fixed negotiation order or a negotiation callback can be defined. For a fixed negotiation order, the agents' doNegotiation
function will be called in the defined order.
clear all;
% should always be called before defining and running another simulation
yalmip('clear');
[hlModel, hlController, T_s_h] = create_hl_ansatz();
hlAgent = Agent('HL', hlModel, hlController, T_s_h, [0 0 0]');
[llModel, llController, T_s_l] = create_ll_ansatz();
llAgent = createAgent('LL', llModel, llController, T_s_l, [0 0 0 0 0]');
sim = Simulation('testsimulation', 24*60);
% set simulation configuration
sim.config.livePlot = true;
% automatically save simulation progress every 10 time steps
sim.config.autoSave = 10;
sim.addAgents(hlAgent, llAgent);
sim.negotiationOrder = {'hl', 'll', 'hl'};
% execute the simulation
sim.runSimulation();
The flow of a simulation is pretty straight forward:
- Check that timesteps of agents are even multiples
- Generate the random number generator intial seed and the seeds for each timestep based on the initial one
- Generate the call order of the agents
- Store seeds
- Simulation loop, step size is maximal step size of the agents
- Set rng seed
- Negotiation loop, if defined
- Call agent's negotiation using
doNegotiation
- Call agent's negotiation using
- Loop over call order of agents (slowest (i.e. longest time step), faster, ..., fastest (i.e. shortest time step))
- Perform agent's simulation step using
doStep
- Perform agent's simulation step using
- Final update of plots
- Store history
Per default, the simulation results will be stored in a sub directory Results/NAME_OF_THE_SIMULATION_YYYYMMDD
of the directory from which the simulation was started, so for the example above Results/testsimulation_20210131
.
Each agent has its own subdirectory, in which the history will be stored (as CSV):
- The state, input and disturbance trajectories (
x.csv
,u.csv
,d.csv
) - The results of each eval function,
eval_NAME.csv
- The results of each cost function,
costs_NAME.csv
- The RNG seeds of the agent
timestep_seeds.csv
- The simulation time vector in agent time
time.csv
- The agent's and controller's configuration in human readable format
config.txt
- The Pareto optimization results, so fronts, parameters etc., in the subdirectory
pareto
, see Pareto Optimization with PARODIS
The virtual history will be stored in the subdirectory virtual
. See the next paragraph for an explanation of the virtual history
Those plots which are configured to be stored will stored in Results/NAME_OF_THE_SIMULATION_YYYYMMDD/plots
.
For more information about the simulation instance, about how results are stored etc, see Simulation.
PARODIS keeps two different histories: The true history and the virtual history.
The true history is the history of actually realised states, inputs and disturbances and thus costs and evaluation results. I.e. that what happens, when the real disturbance measureState
. This history is stored in agent.history
.
Contrary to that, the virtual history, stored in agent.virtualHistory
tracks the history that would have occured if only dPred
was applied in each step, and no state measurement occured.
This is very useful to track differences between predicted and realised behaviour.
Another feature of PARODIS is, that if a simulation should end prematurely (MATLAB crash, solver crash, manual cancellation), you can restore the simulation and continue the simulation from a specified simulation time. This can be done by simply replacing the line
sim.runSimulation();
by
sim.restoreSimulation( directory, k_start );
where directory in this case would be something like "testsimulation_20200907_145030"
, and k_start
the time step from which on to restore the simulation.
Restoring the simulation is done by reading the stored history into the configured objects (entries default to NaN
if they should be missing) and by setting the RNG seeds to the ones of the stored simulation, so that any random number generating functions will behave consistently.
PARODIS supports plotting of states, inputs, disturbances, cost functions and evaluation functions.
If sim.config.livePlot = true
, then plots are updated live during the simulation, while both history as well as the predicted values are plotted.
After the simulation, the figures can be saved as .eps, .png and .fig (default setting) with sim.config.storePlots = 1
.
For these final plots, the prediction horizon is omitted.
Figures can contain multiple subplots and are created as instances of a subclass of the PARODIS-class Figure
.
There is one subclass TimeSeries
shipped with this framework.
When creating the instance (i.e. figure), a name and the number of rows and columns of subplots is defined, e.g.
% define figures
fig1 = TimeSeries("plotMe", 1, 2);
Apperance of the axes are defined by class-functions such as setXLabel()
, see TODO.
Regular plots are defined by the class function addLine()
, which uses the Matlab function stairs()
to plot either states, inputs, disturbances, cost functions or eval functions. Legend entries as well as all valid options can be defined by cell arrays as additional arguments.
Note that by default, x-Axes are linked, such that all subplots share the time limits within the plots.
To prevent linking, use e.g. fig1.setDoLinkAxes(0)
. Furthermore, if no index for a variable is given, then all are taken.
...
% addLine(agent, variable, variableIndex, subplotIndex, legendName, ...
% scenario, optionsReal, optionsPred, leftOrRight)
fig1.addLine(agent1, 'u', [], 1, {'input1', 'input2'}, [], {}, {}, 'left');
fig1.addLine(agent1, 'x', [1], 1, {'x1: blau'}, [], {'Color', 'blue'}, {'Color', 'blue'}, 'right');
fig1.setYLabel(1, 'Power in kW', 'left', {'FontSize', 12, 'Color', 'magenta'})
fig1.setYLabel(1, 'State of Charge', 'right', {'FontSize', 12, 'Color', 'blue'})
% addHeatmap(agent, variable, variableIndexArray, subplotIndex, ...
% scenario, optionsReal, optionsPred)
fig1.addHeatmap(agent1, 'x', [1 3], 2, [], {}, {});
fig1.setYLabel(2, 'EV Number')
% if no index is given, xLabel for all subplots is set
fig1.setXLabel([], 'Time in s', {'FontSize', 12, 'Color', 'blue'})
After all figures are defined, they have to be added to the simulation instance, which can then be started.
...
sim.addPlot(fig1);
% execute the simulation
sim.runSimulation();
- Prerequisites
- Installing PARODIS
- Creating an Ansatz
- Setting up an Agent
- Setting up a Simulation
- Setting up Figures
Pareto Optimization
Modules
- Simulation
- Agent
- Model
- Controllers
- Cost Functions
- Sources
- Plotting
Other