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).






To understand the reason behind the lack of null controllability for the memorytype 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 infintedimensional ODE which, since it has no diffusion term, introduces nonpropagation 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].
Memorytype 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 memorytype 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 nondifferentialbe 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).


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).


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 memorytype equation will be computed throught the minimization problem
Simulations
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:
unzip('https://github.com/DeustoTech/DyConComputationalPlatform/archive/master.zip')
addpath(genpath(fullfile(cd,'DyContoolboxmaster'))
StartDyconPlatform
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
clear;
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(end3),u(end2))*f(1:end2) ; ...
0 ; ...
0 ; ...
f(end1) ; ...
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
end
%  State constraints 
opti.subject_to(Ucas(:,1)==[Y0 ; 0.7 ; 0.7; 1.5 ; 1.5]);
%  Optimization objective 
Cost = (Ucas(1:end4,Nt+1))'*(Ucas(1:end4,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',1e4,'constr_viol_tol',1e3,'compl_inf_tol',1e3);
opti.solver('ipopt',p_opts,s_opts); % set numerical backend
tic
sol = opti.solve(); % actual solve
toc
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$.



Bibliography
[1] U. Biccari and S. Micu, Nullcontrollability properties of the wave equation with a second order memory term, J. Differential Equations, 267(2), 2019, 13761422.
[2] U. Biccari and M. Warma, Nullcontrollability properties of a fractional wave equation with a memory term, Evol. Eq. Control The., 9(2), 2020, 399430.
[3] F. W. ChavesSilva, 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, 500531.