Home

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
Installing PARODIS
Creating an Ansatz
Setting up an Agent
- Evaluation functions
Setting up a Simulation
Setting up Figures

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))](https://render.githubusercontent.com/render/math?math=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)$ .

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)](https://render.githubusercontent.com/render/math?math=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';

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(https://render.githubusercontent.com/render/math?math=[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(https://render.githubusercontent.com/render/math?math=[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

SymbolicController, where the optimization problem is compiled initially as completely parametric YALMIP optimizer 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's optimize() 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 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)](https://render.githubusercontent.com/render/math?math=\Delta x(n|k) = x(n+1|k) - x(n|k)) (or $u$ respectively), again set equally for all scenarios and the entire horizon. 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);

Creating an agent using createAgent, providing a name, an ansatz-handle (which creates a model, a controller and an associated T_s), and a value for the initial state of the system

agent = createAgent(name, @create_ansatz_handle, x0, implicitPrediction);

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.

agent = createAgent('testAgent', @create_ansatz_handle, x0);
% changing configuration parameters
agent.config.solver = 'quadprog';
agent.config.debugMode = true;
agent.config.pareto.calculatePareto = true;
agent.config.pareto.metricFunction = 'RoC';

% 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)
    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');

hlAgent = createAgent('HL', @create_hl_ansatz, [0 0 0]');
llAgent = createAgent('LL', @create_ll_ansatz, [0 0 0 0 0]');

sim = Simulation('testsimulation', 24*60);
% set simulation configuration
sim.config.livePlot = true;
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
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
Final update of plots
Store history

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, T_start );

where directory in this case would be something like "testsimulation_20200907_145030", and T_start can be any time smaller than the time at which the simulation ended.

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'})

% addColorMap(agent, variable, variableIndexArray, subplotIndex, ...
%             scenario, optionsReal, optionsPred)
fig1.addColorMap(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();

Getting Started

Pareto Optimization

Modules

Other

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Home

Getting Started

Overview

Prerequisites

Installing PARODIS

Creating an Ansatz

What is a Scenario

Choosing a Controller

Adding Parameters

Adding Constraints

Cost Functions

Adding Slack Variables independently of Cost Functions

Disturbance Sources

Setting up an Agent

Evaluation functions

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

Setting up Figures

Getting Started

Clone this wiki locally