Work Packages Models involving memory terms and hybrid PDE+ODE systems

Moving control strategy for memory-type equations

Authors: - 20 October 2020

Download Code

This tutorial shows how to implement the moving control strategy for the heat equation in the presence of a memory term.

The standard null controllability problem for the heat equation reads as follows

where $\Omega$ is a bounded domain, $u_0\in L^2(\Omega)$ and the control $f$ acts on the subdomain $\omega\subset\Omega$. The objective is to steer the system from $u_0$ to zero in time $T$,

It is known that the null controllability of the heat equation holds in any positive time $T>0$ and for any open control region $\omega$. Nevertheless, the situation changes when in the model appears a memory term:

As a matter of fact, when the equation involes a memory term, it is impossible to achive null controllability in the classical sense, unless the control region coincides with the entire domain ($\omega=\Omega$).

This fact can be easily seen when applying a LQR control on both systems (\ref{heatinterior}) and (\ref{heatinterior_memory}). In the former case (without memory) the minimum of the LQR functional is an optimal control stabilizing the dynamics (Video 1), whereas when the memory enters into the model this minimum leaves the system far from the equilibrium (Vdeo 2).

Video 1: LQR in the heat equation (Equation \ref{heatinterior})
Video 2: LQR in the heat equation with memory (Equation \ref{heatinterior_memory})
We can see the time evolution of the two models with a control region, $\chi_{\omega}$, located in the position of the blue ellipsoid in the simulation. We observe that the LQR control is effective for the heat equation but not for the heat equation with memory.

To understand the reason behind the lack of null controllability for the memory-type equation (\ref{heatinterior_memory}), we can introduce the new variable

In this way, (\ref{heatinterior_memory}) is transformed into a coupled PDE/ODE system:

In (\ref{heatinterior_memory_z}), the equation $z_t = u$ is an infinte-dimensional ODE which, since it has no diffusion term, introduces non-propagation phenomena in the system. The non propagating components corresponding to this ODE are unable to reach the control region $\omega$ and, therefore, cannot be controlled.

A natural remedy to this issue is to operate with a moving control strategy. If the ODE components of the system cannot reach the control region $\omega$, then is the control region who needs to reach these components. Hence, our control region would be function of the time variable, $\omega = \omega(t)$, and the corresponding model will be

This moving control strategy has been proposed and succesfully applied in several contributions of our team [1,2,3,4].

Memory-type equations are one of the core topics of WP5 of the DyCon ERC project. An abridged presentation of this working package can be found at the following link.

Numerical implementation of the moving control strategy

We present here the numerical implementation of the moving control strategy for the memory-type heat equation (\ref{heatinterior_memory_z_moving}). All the simulations will be in $2D$.

STEP 1: construction of the control region

The first step is to construct the moving control region $\omega(t)$, which requires to design the function $\chi_{\omega(t)}$. To this end, we made the assumption that the size and orientation of $\omega(t)$ do not change during the time interval $[0,T]$. Hence, once fixed the size and shape of the control region, we will obtain a function depending only on a point $x\in \Omega$.

In this simulation, we create the control region as a square which moves in the domain $\Omega$. This square can be constructed with a function $W(x)$ defined as

where $H(x)$ is the Heaviside function. Nevertheless, the function $W$ constructed in this way would be non-differentialbe at the points $x=a$ and $x=b$. This would represent a problem when implementing the gradient method to compute the control. To bypass this issue, in the definition of $W(x,a,b)$ we will use a smooth approximation of $H$ (See Figure 1).

Figure 1: Graphicla representation of $W(x,a,b)$

Moreover, the generalization of this function to dimension two is immediate

In this way, given $(x_{min}, x_{max}, y_{min}, y_{max})$, we obtain the corresponding characteristic function which only depends on the position of a single point.

STEP 2: construction of the dynamics

Since $\chi_{\omega}$ can move and, moreover, we want to obtain the optimal trajectory for the control, we will incorporate in our system a moving particle whose position and velocity are denoted by $\textbf{d} = (d_x,d_y)$ and $\textbf{v} = (v_x,v_y)$, respectively (Video 3).

Video 3: Representación gráfica de $W_{2D}$

Besides, we will add a control $\textbf{g}(t) = (g_x(t),g_y(t))$ which will act as an exterior force moving the subset $\chi_{\omega}$.

In this way, our model takes the final form

STEP 3: optimal control problem

Following a standard optimal control strategy, the control function for our memory-type equation will be computed throught the minimization problem


All the simulations in this tutorial have been performed using the DyCon Toolbox and CasADi. Since CasADi is incorporated into the DyCon toolbox, we can install it in the following way:


We shall now write the discrete version of equation (\ref{heatinterior_memory_z_all}):

Secondly, let us create the mesh variable

% Create the mesh variables
Ns = 7;
Nt = 15;
xline = linspace(-1,1,Ns);
yline = linspace(-1,1,Ns);
[xms,yms] = meshgrid(xline,yline);

and the matrix $A$ containing the $2D$ Laplacian and the ODE dynamics

A  = FDLaplacial2D(xline,yline);

Atotal = zeros(2*Ns^2+4,2*Ns^2+4);
Atotal( 1:Ns^2  , 1:Ns^2 ) = A;
Atotal( Ns^2+1 : 2*Ns^2   ,   1    :  Ns^2   )  =  eye(Ns^2);
Atotal(    1   :  Ns^2    , Ns^2+1 : 2*Ns^2  )  =  50*eye(Ns^2); % z = 50*y

RumbaMatrixDynamics = [0 0 1 0; ...
                       0 0 0 1; ...
                       0 0 0 0; ...
                       0 0 0 0 ];

Atotal(2*Ns^2+1:end,2*Ns^2+1:end) = RumbaMatrixDynamics;
Atotal = sparse(Atotal);

Finally, we create the function $B(\textbf{d}) = B(x_d,y_d)$ for the dynamics of the moving control

% We create the B() function
xwidth = 0.3;
ywidth = 0.3;
B = @(xms,yms,xs,ys) WinWP05(xms,xs,xwidth).*WinWP05(yms,ys,ywidth);
Bmatrix =  @(xs,ys) [diag(reshape(B(xms,yms,xs,ys),1,Ns^2)) ;zeros(Ns^2)];

We now have everything we need to construct and solve the optimal control problem

opti = casadi.Opti();  % CasADi optimization structure

% ---- Input variables ---------
Ucas = opti.variable(2*Ns^2+4,Nt+1); % state trajectory
Fcas = opti.variable(Ns^2+2,Nt+1);   % control

% ---- Dynamic constraints --------
dUdt = @(y,f) Atotal*u+ [Bmatrix(u(end-3),u(end-2))*f(1:end-2) ; ...
                                         0                     ; ...
                                         0                     ; ...
                                     f(end-1)                  ; ...
                                     f(end)                    ]; 

% -----Euler backward method-------
for k=1:Nt % loop over control intervals
   y_next = Ucas(:,k) + (T/Nt)*dUdt(Ucas(:,k+1),dUdt(:,k+1)); 
   opti.subject_to(Ucas(:,k+1)==y_next); % close the gaps

% ---- State constraints --------
opti.subject_to(Ucas(:,1)==[Y0 ; 0.7 ; 0.7; -1.5 ; -1.5]);

% ---- Optimization objective  ----------
Cost = (Ucas(1:end-4,Nt+1))'*(Ucas(1:end-4,Nt+1));

opti.minimize(Cost); % minimizing L2 at the final time

% ---- initial guesses for solver ---
opti.set_initial(Ucas, Unum_free);
opti.set_initial(Fcas, 0);

% ---- solve NLP              ------
p_opts = struct('expand',false);
s_opts = struct('acceptable_tol',1e-4,'constr_viol_tol',1e-3,'compl_inf_tol',1e-3);
opti.solver('ipopt',p_opts,s_opts); % set numerical backend
sol = opti.solve();   % actual solve

The following video shows the results of our simulations, which allowed to compute an effective moving control steering the dynamics to rest at time $T$.

Video 3


[1] U. Biccari and S. Micu, Null-controllability properties of the wave equation with a second order memory term, J. Differential Equations, 267(2), 2019, 1376-1422.

[2] U. Biccari and M. Warma, Null-controllability properties of a fractional wave equation with a memory term, Evol. Eq. Control The., 9(2), 2020, 399-430.

[3] F. W. Chaves-Silva, X. Zhang and E. Zuazua, Controllability of evolution equations with memory, SIAM J. Control Optim., 55(4), 2017, 2437–2459.

[4] Q. Lü, X. Zhang and E. Zuazua, Null controllability for wave equations with memory, J. Math. Pures Appl., 108 (4), 2017, 500-531.