Work Packages Control of parameter dependent problems

Simultaneous Control with DyCon Toolbox

Author: - 01 April 2020

Download Code

In this tutorial we will present a simultaneous control problem in a linear system dependent on parameters. We will use the MATLAb DyCon Toolbox library.

Setup

In this tutorial we need DyCon Toolbox, to install it we will have to write the following in our MATLAB console:

unzip('https://github.com/DeustoTech/DyCon-Computational-Platform/archive/master.zip')
addpath(genpath(fullfile(cd,'DyCon-toolbox-master')))

Problem definition

The simultaneous control problem is defined as:

subject to:

where:

Numerical Implementation

First we import the CasAdi library in order to create symbolic variables

import casadi.*
%% In this case $\nu_i$ are
M = 50;
nu = linspace(1,6,M);

Then we create the matrices of linear dynamics.

[A,B] = GenMatSim(nu);

With these matrices we create the ode object

Nt = 500;T  = 0.8;
tspan = linspace(0,T,Nt);
%% create linear dynamic
iode = linearode(A,B,tspan);
%% set initial condition
Y0 = ones(2, 1);
iode.InitialCondition = repmat(Y0,M,1);

then we create the optimal control problem

% Get Symbolical variable
Ys  = iode.State.sym;
Us  = iode.Control.sym;
ts  = SX.sym('t'); %% <= Create a symbolical time
% Set Target
YT = zeros(2, 1);YT = repmat(YT,M,1);
%
PathCost  = casadi.Function('L'  ,{ts,Ys,Us},{ (1/2)*(Us'*Us)           });
FinalCost = casadi.Function('Psi',{Ys}      ,{  1e7*((Ys-YT).'*(Ys-YT)) });
% Create the optimal control
iocp = ocp(iode,PathCost,FinalCost);

Solve Optimal Control Problem with ipopt solver

U0 = ZerosControl(iode);
[Uopt ,Yopt] =  IpoptSolver(iocp,U0);
This is Ipopt version 3.12.3, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:   199700
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      600

Total number of variables............................:    50500
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:    50000
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.0000000e+09 1.12e-02 7.88e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  4.0934345e+06 1.26e-13 5.63e-09  -1.0 1.66e+02    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:   2.0467172272726700e+01    4.0934344545453396e+06
Dual infeasibility......:   5.6305735629536002e-09    1.1261147125907198e-03
Constraint violation....:   1.2612133559741778e-13    1.2612133559741778e-13
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   5.6305735629536002e-09    1.1261147125907198e-03


Number of objective function evaluations             = 2
Number of objective gradient evaluations             = 2
Number of equality constraint evaluations            = 2
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 1
Total CPU secs in IPOPT (w/o function evaluations)   =      1.391
Total CPU secs in NLP function evaluations           =      0.027

EXIT: Optimal Solution Found.
               t_proc [s]   t_wall [s]    n_eval
       nlp_f     0.000381      0.00038         2
       nlp_g       0.0042       0.0042         2
  nlp_grad_f      0.00167      0.00167         3
  nlp_hess_l      0.00331      0.00332         1
   nlp_jac_g        0.026        0.026         3
      solver         1.48         1.36         1
Elapsed time is 3.853462 seconds.

Compute Free solution

Yfree = solve(iode,U0*0);
Yfree = full(Yfree);

Visualization

fig = figure;
fig.Units = 'norm';fig.Position = [0.1 0.1 0.6 0.5];
%%
plotSimu(tspan,Yfree,Yopt,Uopt,M)

Figure 1. The different colors represent the dynamic system under different parameters. It can be seen how the same control is obtained acting for all the dynamic systems is capable of driving the systems to the target.