Skip to content
jengstud edited this page Aug 6, 2021 · 19 revisions

Getting Started

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 $x$, a control input $u$ and a disturbance input $d$. The model contains symbolic expressions of state trajectory for each scenario x, symbolic expression for each disturbance d and symbolic expression for the control input u. A model is defined by a function that returns an ordinary difference equation for each time step within the agent's prediction horizon.

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.


Overview

Prerequisites

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 PARODIS

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

Creating an Ansatz

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 $x_s(n+1|k) = f_n(x_s(n|k), u(n|k), d_s(n|k))$, to build the symbolic expressions for the state trajectory in multiple scenarios. Since PARODIS defines a prediction horizon by the time steps within the horizon, each time step has it's own associated ODE $f_n(\cdot)$.

If the argument implicitPrediction is set to false, the state trajectory $x_s(n|k)$ will be expressed directly in terms of $x_s(n+1|k) = f_n(x_s(n|k), u(n|k), d_s(n|k))$, so it will be a symbolic expression of $u(n|k)$ and $d_s(n|k)$. If it is set to true, the prediction will be expressed implicitly, meaning that the prediction dynamics will be represented by equality constraints $x_s(n+1|k) == f_n(x_s(n|k), u(n|k), d_s(n|k))$. Depending on your model, it may make sense to use the implicit prediction form, as it can help with problem conditioning.

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 $x(k+1) = Ax(k)+Bu(k)+Sd(k)$, you can let your model function additionally return the matrices $A$, $B$, and $S$. PARODIS will use these to speed up internal calculations:

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.

What is a Scenario

An important expression in PARODIS language is the scenario. A scenario $s$ is the combination of one set of disturbance predictions $[d_s(0|k),\ \dots,\ d_s(N_{pred}-1|k)]$ that is generated at each time step $k$, the resulting state trajectory $[x_s(0|k),\ \dots,\ x_s(N_{pred}|k)]$, and a set of parameter values. In PARODIS, it is assumed that the same optimal input $u$ is applied in each 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.

Choosing a Controller

PARODIS offers three different controllers

  1. SymbolicController, where the optimization problem is compiled initially as completely parametric YALMIP optimizer object, i.e. the problem remains fully parametrized at all times
  2. ExplicitController where the optimization problem is recreated at each time step, parameters are always directly replaced by their values and YALMIP's optimize() function is used
  3. 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 the ExplicitController and then creates a YALMIP optimizer 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.

Adding Parameters

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.

Adding Constraints

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 for du(0|k) = u(0|k) - u(k-1) the input u(k-1) of the previous time step is automatically provided by PARODIS. At simulation start k = 0, PARODIS will either assume u(k-1) to be zero, or use an initial value provided to the Agent constructor (see Agent). The argument dT 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

Cost Functions

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

Adding Slack Variables independently of Cost Functions

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.

Disturbance Sources

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.


Setting up an Agent

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

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. If addEvalFunction( 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).


Setting up a Simulation

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:

  1. Check that timesteps of agents are even multiples
  2. Generate the random number generator intial seed and the seeds for each timestep based on the initial one
  3. Generate the call order of the agents
  4. Store seeds
  5. Simulation loop, step size is maximal step size of the agents
    1. Set rng seed
    2. Negotiation loop, if defined
      1. Call agent's negotiation using doNegotiation
    3. Loop over call order of agents (slowest (i.e. longest time step), faster, ..., fastest (i.e. shortest time step))
      1. Perform agent's simulation step using doStep
  6. Final update of plots
  7. Store history

Simulation Results

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.

History and virtual history

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 $d$ is applied to the system and/or when the history was manipulated in 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.

Restoring a Simulation

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.

Setting up Figures

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();
Clone this wiki locally