# Control for a semilinear heat equation and analogies with a collective behavior model

Author: - 15 March 2020

Download Code

In this tutorial we will apply the DyCon toolbox to find a control to the semi-discrete semi-linear heat equation.

where $N^2A$ is the discretization of the Laplacian in 1d in $N$ nodes. We are looking for a control that after time $T$ steers the system near zero.

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


### Methodology

In order to do so we will frame the problem as a minimization of a functional and we will apply gradient descent to find it. Note that the convexity of the functional is not proven, therefore, we will obtain a local minima for the functional. The functional considered will be:

where by $| \cdot |_{L^2}$ we understand the discrete $L^2$ norm.

Once this control is computed for a certain N, we will think on a dynamical system that models an opinion dynamics with $N$ agents communicating through a chain.

\begin{equation}\label{m1}y_t-\frac{1}{N}Ay=G(y)+Bv\end{equation}

The goal will be to compute also the control $v$ thinking model \eqref{m1} as if it was a semidiscretization of a heat equation with diffusivity $\frac{1}{N^3}$. Furthermore we will also compute the control for model \eqref{m1} with the non-linearity being non-homogeneous on $N$ and a time horizon being $T_N=N^3T$.

clear all
import casadi.*


Definition of the time

ts = SX.sym('t');
%% Discretization of the space
N = 30;
xi = 0; xf = 1;
xline = linspace(xi,xf,N);


Interior Control region between 0.5 and 0.8

w1 = 0.5; w2 = 0.8;
D = 1;

[A,B,count] = GetABmatrix(N,w1,w2,D);


we define symbolically the vectors of the state and the control

Ys = SX.sym('y',[N 1]);
Us = SX.sym('u',[count 1]);


We create the ODE object Our ODE object will have the semi-discretization of the semilinear heat equation. We set also initial conditions, define the non linearity and the interaction of the control to the dynamics.

Initial condition

Y0 = 2*sin(pi*xline');


Diffusion part: the discretization of the 1d Laplacian

Definition of the non-linearity $\partial_y[-5\exp(-y^2)]$

G = casadi.Function('NLT',{Ys},{10*Ys.*exp(-Ys.^2)});%%


and we define the part of the dynamics corresponding to the nonlinearity

Putting all the things together

F = casadi.Function('F',{ts,Ys,Us},{A*Ys + B*Us + G(Ys)});
%% Creation of the ODE object
%% Time horizon
T = 1;


We create the ODE-object and we change the resolution to $dt=0.01$ in order to see the variation in a small time scale. We will get the values of the solution in steps of size odeEqn.dt, if we do not care about modifying this parameter in the object, we might get the solution in certain time steps that will hide part of the dynamics.

Nt = 100;
tspan = linspace(0,T,Nt);
ipde = semilinearpde1d(Ys,Us,A,B,G,tspan,xline);
ipde.InitialCondition = Y0;


We solve the equation and we plot the free solution applying solve to odeEqn and we plot the free solution.

Yfree = solve(ipde,ZerosControl(ipde));
Yfree = full(Yfree);

figure;
surf(Yfree','EdgeColor','none');
title('Free Dynamics')
xlabel('space discretization')
ylabel('Time')
yticks([1 Nt])
yticklabels([tspan(1) tspan(end)]) We create the object that collects the formulation of an optimal control problem by means of the object that describes the dynamics odeEqn, the functional to minimize Jfun and the time horizon T

L   = Function('L'  ,{ts,Ys,Us},{ Us.'*Us  });
Psi = Function('Psi',{Ys}      ,{ 1e6*(Ys.'*Ys) });

iocp = ocp(ipde,L,Psi);


We apply the steepest descent method to obtain a local minimum (our functional might not be convex).

U0 =ZerosControl(ipde);
[OptControl ,OptState]  = 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...:    19236
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     3900

Total number of variables............................:     3900
variables with only lower bounds:        0
variables with lower and upper bounds:        0
variables with only upper bounds:        0
Total number of equality constraints.................:     3000
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  3.0000000e+07 9.05e+00 9.67e-01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
1  2.4544321e+05 1.22e-01 7.92e+00  -1.0 3.10e+01  -2.0 1.00e+00 1.00e+00f  1
2  1.8219937e+06 4.24e-03 1.73e+00  -1.0 1.31e+01  -1.6 1.00e+00 1.00e+00h  1
3  7.1220860e+04 1.67e-02 2.92e+00  -1.0 9.12e+00  -2.1 1.00e+00 1.00e+00f  1
4  8.9356500e+03 3.97e-04 8.51e-02  -1.0 2.41e+00  -2.5 1.00e+00 1.00e+00f  1
5  3.1075917e+03 1.06e-04 3.19e-03  -2.5 2.96e+00  -3.0 1.00e+00 1.00e+00f  1
6  1.3016546e+03 5.20e-05 1.19e-03  -3.8 3.60e+00  -3.5 1.00e+00 1.00e+00h  1
7  4.5871538e+02 5.64e-05 4.84e-04  -5.7 4.41e+00  -4.0 1.00e+00 1.00e+00h  1
8  1.6346578e+02 6.55e-05 1.92e-04  -5.7 5.26e+00  -4.4 1.00e+00 1.00e+00h  1
9  9.7049320e+01 1.20e-04 3.78e-05  -5.7 3.10e+00  -4.9 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
10  4.0069187e+01 9.26e-03 1.98e-05  -5.7 1.48e+01    -  1.00e+00 5.00e-01h  2
11  1.6748754e+01 4.28e-04 2.92e-06  -5.7 4.34e+00    -  1.00e+00 1.00e+00h  1
12  1.6371256e+01 4.58e-06 9.87e-09  -5.7 1.88e-01    -  1.00e+00 1.00e+00h  1
13  1.6367177e+01 1.03e-11 1.11e-13  -8.6 7.84e-04    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 13

(scaled)                 (unscaled)
Objective...............:   8.1835886066017175e-04    1.6367177213203433e+01
Dual infeasibility......:   1.1096830919432588e-13    2.2193661838865175e-09
Constraint violation....:   1.0345974077452524e-11    1.0345974077452524e-11
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   1.0345974077452524e-11    2.2193661838865175e-09

Number of objective function evaluations             = 19
Number of objective gradient evaluations             = 14
Number of equality constraint evaluations            = 19
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 14
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 13
Total CPU secs in IPOPT (w/o function evaluations)   =      0.436
Total CPU secs in NLP function evaluations           =      0.057

EXIT: Optimal Solution Found.
t_proc [s]   t_wall [s]    n_eval
nlp_f     0.000826     0.000828        19
nlp_g      0.00916      0.00917        19
nlp_grad_f      0.00139      0.00139        15
nlp_hess_l       0.0175       0.0174        13
nlp_jac_g       0.0333       0.0333        15
solver        0.507        0.501         1
Elapsed time is 0.636870 seconds.


plotSHE(OptState,OptControl,Yfree,YT,xline,tspan,Nt)  ### Collective behavior dynamics

Now we apply the same procedure for the collective behavior dynamics. We will employ a function that does the algorithm explained before for the semilinear heat equation having the chance to set a diffusivity constant.

For the simulation of the model in collective behavior we will employ a diffusivity $D=\frac{1}{N^3}$.

D = 1/(N^3);
[A,B,count] = GetABmatrix(N,w1,w2,D);
%%
G = casadi.Function('NLT',{Ys},{10*Ys.*exp(-Ys.^2)});%%
%%
T = 1;  Nt = 100;
tspan = linspace(0,T,Nt);
%%
ipde = semilinearpde1d(Ys,Us,A,B,G,tspan,xline);
ipde.InitialCondition = Y0;
%%
iocp = ocp(ipde,L,Psi);  U0 =ZerosControl(ipde);
%%[OptControl ,OptState]  = ArmijoGradient(iocp,U0,'MaxIter',200);
[OptControl ,OptState]  = 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...:    19236
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     3900

Total number of variables............................:     3900
variables with only lower bounds:        0
variables with lower and upper bounds:        0
variables with only upper bounds:        0
Total number of equality constraints.................:     3000
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  3.0000000e+07 1.00e+00 9.72e-01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
1  1.1211823e+08 6.93e-02 1.07e+01  -1.7 1.75e+00   0.0 1.00e+00 1.00e+00h  1
2  1.0201471e+08 2.53e-02 2.23e+00  -1.7 3.50e+00  -0.5 1.00e+00 1.00e+00f  1
3  8.8642262e+07 4.18e-02 5.31e-01  -1.7 2.13e+02    -  1.00e+00 1.00e+00f  1
4  8.2505158e+07 1.39e+00 3.13e+00  -1.7 2.72e+03    -  1.00e+00 1.00e+00F  1
5  7.9212105e+07 9.68e+00 5.20e+00  -1.7 1.38e+03    -  1.00e+00 1.00e+00f  1
6  7.8291279e+07 6.94e-02 2.84e+00  -1.7 2.05e+03    -  1.00e+00 1.00e+00f  1
7  7.8090072e+07 1.57e+00 4.98e-01  -1.7 7.87e+02    -  1.00e+00 1.00e+00f  1
8  7.8026943e+07 1.54e-01 8.25e-02  -1.7 3.33e+02    -  1.00e+00 1.00e+00f  1
9  7.8491521e+07 4.41e-02 7.64e+00  -2.5 2.10e+00  -1.0 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
...
...
54  7.8022238e+07 4.03e-07 8.33e-10  -5.7 2.15e-02    -  1.00e+00 1.00e+00h  1
55  7.8022238e+07 9.83e-12 1.82e-12  -8.6 1.02e-04    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 55

(scaled)                 (unscaled)
Objective...............:   3.9011119130960064e+03    7.8022238261920124e+07
Dual infeasibility......:   1.8189894035458565e-12    3.6379788070917130e-08
Constraint violation....:   9.8296371042749797e-12    9.8296371042749797e-12
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   9.8296371042749797e-12    3.6379788070917130e-08

Number of objective function evaluations             = 111
Number of objective gradient evaluations             = 56
Number of equality constraint evaluations            = 111
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 56
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 55
Total CPU secs in IPOPT (w/o function evaluations)   =      1.374
Total CPU secs in NLP function evaluations           =      0.308

EXIT: Optimal Solution Found.
t_proc [s]   t_wall [s]    n_eval
nlp_f      0.00624      0.00627       111
nlp_g       0.0634       0.0633       111
nlp_grad_f      0.00762      0.00768        57
nlp_hess_l       0.0868       0.0846        55
nlp_jac_g        0.149        0.149        57
solver         1.74         1.73         1
Elapsed time is 1.874814 seconds.


plotSHE(OptState,OptControl,Yfree,YT,xline,tspan,Nt)  Now we will change also the time horizon and we will incorporate a non-homogeneous non-linearity, we will just divide the non-linearity $G$ by $N^3$

D = 1/(N^2);
[A,B,count] = GetABmatrix(N,w1,w2,D);
%%
G = casadi.Function('NLT',{Ys},{10*Ys.*exp(-Ys.^2)/N^2} );%%
%%
T = N^2;  Nt = 100;
tspan = linspace(0,T,Nt);
%%
ipde = semilinearpde1d(Ys,Us,A,B,G,tspan,xline);
ipde.InitialCondition = Y0;
%%
iocp = ocp(ipde,L,Psi);

U0 =ZerosControl(ipde);
[OptControl ,OptState]  = 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...:    19236
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     3900

Total number of variables............................:     3900
variables with only lower bounds:        0
variables with lower and upper bounds:        0
variables with only upper bounds:        0
Total number of equality constraints.................:     3000
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  3.0000000e+07 9.05e+00 1.17e+02  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
1  2.8754020e+04 1.48e-01 4.01e+00  -1.0 2.73e+00   0.0 1.00e+00 1.00e+00f  1
2  3.0839145e+03 1.91e-02 3.08e-01  -1.0 9.19e-01  -0.5 1.00e+00 1.00e+00h  1
3  9.2703794e+02 1.17e-03 2.31e-02  -1.7 2.08e-01  -1.0 1.00e+00 1.00e+00h  1
4  4.2707593e+02 2.49e-03 9.65e-03  -3.8 2.61e-01  -1.4 1.00e+00 1.00e+00h  1
5  2.2851845e+02 1.36e-03 2.60e-03  -3.8 2.10e-01  -1.9 1.00e+00 1.00e+00h  1
6  1.5158663e+02 1.25e-03 9.80e-04  -3.8 2.38e-01  -2.4 1.00e+00 1.00e+00h  1
7  1.0229631e+02 2.37e-03 4.26e-04  -5.7 3.11e-01  -2.9 1.00e+00 1.00e+00h  1
8  8.6220145e+01 1.84e-01 3.73e-04  -5.7 2.84e+01    -  1.00e+00 1.25e-01h  4
9  1.6663658e+00 6.98e-02 4.71e-05  -5.7 2.20e+00    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
10  1.1303218e-01 1.61e-02 9.28e-07  -5.7 4.53e-01    -  1.00e+00 1.00e+00h  1
11  3.7959122e-02 7.79e-03 2.54e-08  -5.7 2.84e-01    -  1.00e+00 1.00e+00h  1
12  2.1316195e-02 1.77e-03 4.98e-09  -5.7 1.39e-01    -  1.00e+00 1.00e+00h  1
13  1.8447375e-02 6.30e-06 5.26e-11  -5.7 8.15e-03    -  1.00e+00 1.00e+00h  1
14  1.8440851e-02 2.26e-09 9.54e-15  -8.6 1.61e-04    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 14

(scaled)                 (unscaled)
Objective...............:   9.2204254694796125e-07    1.8440850938959225e-02
Dual infeasibility......:   9.5381182640930876e-15    1.9076236528186175e-10
Constraint violation....:   2.2559639156760625e-09    2.2559639156760625e-09
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   2.2559639156760625e-09    2.2559639156760625e-09

Number of objective function evaluations             = 19
Number of objective gradient evaluations             = 15
Number of equality constraint evaluations            = 19
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 15
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 14
Total CPU secs in IPOPT (w/o function evaluations)   =      0.415
Total CPU secs in NLP function evaluations           =      0.065

EXIT: Optimal Solution Found.
t_proc [s]   t_wall [s]    n_eval
nlp_f     0.000898     0.000903        19
nlp_g      0.00904      0.00905        19
nlp_grad_f       0.0017      0.00169        16
nlp_hess_l       0.0214       0.0214        14
nlp_jac_g       0.0354        0.035        16
solver        0.489        0.487         1
Elapsed time is 0.616536 seconds.


plotSHE(OptState,OptControl,Yfree,YT,xline,tspan,Nt)  