#### 1 Introduction

In the post IpOpt and AMPL use to solve time optimal control problems, we explain how to use *IpOpt* and *AMPL* in order to solve control problems with control constraints and possibly some state constraints.

In the present post, we are going to present a numerical development in order to find the minimal controllability time for the discretized heat equation with unilateral (non-negative) control constraint. More precisely, we consider the controllability of a discretized version of the 1D heat equation under nonnegative Dirichlet control constraint. The infinite dimensional system, we consider is

where $\psi$ is the state to control $\psi^0$ the initial state and $u$ the Dirichlet control. The continuous version of this problem has been analysed in [2] and it has been shown that for every initial state $\psi^0\in L^2(0,1)$ and every positive constant target $\psi^1$, the controllability of this system under the control constraint \begin{equation}\label{eq:ContConst} u(t)\geqslant0 (t\geqslant 0), \end{equation} requires a positive waiting time as soon as the target state $\psi^1$ is different from the initial state $\psi^0$.

In this post, we consider a spatially discretized version with has the form

where the matrices $A$ and $B$ are

where $n+1>2$ is the number of discretization points and $[y]_i(t)$ (the $i^\text{th}$ component of $y(t)\in\mathbb{R}^n$) stands for $\psi(t,(i-1)/n)$. To obtain the above finite dimensional system, we have used centred finite differences.

The aim is then to minimize the time $T$ such that there exist a control $u:[0,T]\to\mathbb{R}_+$ steering the solution of (\ref{eq:SysGen}) from $y^0\in\mathbb{R}^n$ to $y^1\in\mathbb{R}^n$ in time $T$. For the sake of simplicity, we assume that $y^0$ and $y^1$ are constant vectors of $\mathbb{R}^n$, that is to say that

for some $u^0\in\mathbb{R}$ and $u^1\in\mathbb{R}$ Since the control $u$ shall satisfy the constraint (\ref{eq:ContConst}), it is easy to see, using a discrete version of the comparison principle, that in order to have the existence of a solution, then we need $u^1>0$. Note that, using again the comparison principle, we can also show that if $u^0\geqslant 0$ then, whatever the control $u\geqslant 0$ is, the solution of (\ref{eq:SysGen}) always satisfies $y(t)\geqslant 0$.

To conclude, the optimization problem we aim to solve is

where $y^0$ and $y^1$ are given by (\ref{eq:y0y1}) for some $u^0\in\mathbb{R}$ and $u^1\in\mathbb{R}_+^*$, and where $y(t;u;y^0)$ stands for the solution of (\ref{eq:SysGen}) at time $t$, with control $u$ and initial condition $y^0$.

In this post, in order to find numerically the minimal controllability time and a time-optimal control, we are going to use the abstract result presented in the next paragraph. Another way could be to introduce, in addition to the control constraint $u(t)\geqslant 0$, the additional control constraint $u(t)\leqslant M$, and let $M$ increase to $+\infty$. Such an approach has been introduced in the previous post IpOpt and AMPL use to solve time optimal control problems.

#### 2 Abstract result

It is shown in [1] that at the minimal time $\underline{T}(y^0,y^1)$, defined by (\ref{eq:optim0}), there exists a non-negative control $u$ in the class of Radon measure. Furthermore, it is proved in this article that at the minimal controllability time, there exists one and only one non-negative Radon measure and this time optimal control is a convex sum of at most $\lfloor (n+1)/2 \rfloor$ Dirac masses. In other words, the non-negative time optimal control takes the form:

with $N=\lfloor (n+1)/2 \rfloor$, $0\leqslant \alpha_k$ and $0\leqslant t_1\leqslant\dots\leqslant t_N\leqslant T$. For a control of this form, the solution of (\ref{eq:SysGen}) at time $T$, with initial condition $y^0$, is:

Consequently, the minimal time $\underline{T}(y^0,y^1)$ given by (\ref{eq:optim0}) is a minimizer of:

#### 3 Numerical implementation

The optimization problem (\ref{eq:optim1}) is not easy to solve directly, mainly because of the presence of a matrix exponential. Thus, instead we will solve the discretized heat equation on each time interval $(t_k,t_{k+1})$, with $t_0=0$ and $t_{N+1}=T$. Let us then set $\tau_k=t_{k+1}-t_k$ for every $k\in{0,\dots,N}$. We have $T=\sum_{k=0}^N\tau_k$ and the optimization problem (\ref{eq:optim1}) also writes

In the above constraints, $y_k(\tau)=y(\tau+t_k)$ for every $\tau\in(0,\tau_k)$, where $y$ is the solution of (\ref{eq:SysGen}), with control $u=\sum_{k=1}^N\alpha_k\delta_{t_k}$ and initial condition $y(0)=y^0$. Notice that since $\alpha_k$ is only constrained to be non-negative, and since the vector $B$ (given by (\ref{eq:ABheat})) is of the form $(0,\dots,0,b_n)^\top$, with $b_n\geqslant 0$, the constraint $y_{k+1}(0)=y_k(\tau_k)+\alpha_{k+1}B$ can be expressed as

Consequently, the parameters $\alpha_1,\dots,\alpha_N$ can be forgotten.

In order to numerically compute $y_k(t)=e^{tA}y_k(0)$ for $t\in[0,\tau_k]$, we are going to use the Crank-Nicholson method. More precisely, given $N_t\in\mathbb{N}^*$, we approximate $y_k(j\tau_k/N_t)$ by $y_k^j$, with $y_k^j$ solution of

All in all, the fully discretized optimization problem, is

subject to the constraints

Note that instead of imposing $[y_{N}^{N_t}]_i=u^1$, we chose to relax this condition in the constraint (\ref{eq:optD6}). This choice has been a maid since we do not expect that the numerical solution exactly reaches the target $y^1$. In practice, we will take, $\varepsilon=1/(n\, N_t)$.

In term of *AMPL* language, the above constrained optimization problem is

# define parametersparamn = 20; # number of spatial discretization pointsparamNt = 400; # number of time discretization pointsparamN = floor((n+1)/2); # maximal number of control impulsesparameps= 1/(n*Nt); # relaxation parameter $\varepsilon}$ in (\ref{eq:optD6})paramu0 = 1; # parameters $u^0\text{ and }u^1}$, see (\ref{eq:y0y1})paramu1 = 5; # define variablesvary {k in 0..N, j in 0..Nt, i in 1..n}; # unknowns $[y_k^j]_i$vartau {k in 0..N} >=0; # unknown impulse times $\tau_k}$ subject to (\ref{eq:optD2}) # objective function, see (\ref{eq:optD1})minimizeT: sum {k in 0..N} tau[k]; # define the constraints # constraint (\ref{eq:optD7}) for interior pointssubject toy_dyn {k in 0..N, j in 1..Nt, i in 2..n-1}: y[k,j,i]-y[k,j-1,i] = (n+1)^2*(y[k,j,i-1]-2*y[k,j,i]+y[k,j,i+1] +y[k,j-1,i-1]-2*y[k,j-1,i]+y[k,j-1,i+1])/2*tau[k]/Nt; # constraint (\ref{eq:optD7}) for i=nsubject toy_dyn_Dirichlet {k in 0..N, j in 1..Nt}: y[k,j,n]-y[k,j-1,n] = (n+1)^2*(y[k,j,n-1]-2*y[k,j,n]+y[k,j-1,n-1]-2*y[k,j-1,n])/2*tau[k]/Nt; # constraint (\ref{eq:optD7}) for i=1subject toy_dyn_Neuman {k in 0..N, j in 1..Nt}: y[k,j,1]-y[k,j-1,1] = (n+1)^2*(y[k,j,2]-y[k,j,1]+y[k,j-1,2]-y[k,j-1,1])/2*tau[k]/Nt; # initial condition, constraint (\ref{eq:optD3})subject toy_init {i in 1..n}: y[0,0,i] = u0; # terminal condition, constraint (\ref{eq:optD6})subject toy_end1 {i in 1..n}: y[N,Nt,i] >= u1-eps;subject toy_end2 {i in 1..n}: y[N,Nt,i] <= u1+eps; # constraint (\ref{eq:optD4})subject tocontinuity {k in 0..N-1, i in 1..n-1}: y[k+1,0,i]=y[k,Nt,i]; # constraint (\ref{eq:optD5})subject tojump {k in 0..N-1}: y[k+1,0,n]>=y[k,Nt,n]; # solve the problem with IpOptoptionsolver ipopt;optionipopt_options "max_iter=100000 linear_solver=mumps halt_on_ampl_error yes";solve; # display parameters and solutionprintf: "u0 = %24.16e\n", u0;printf: "u1 = %24.16e\n", u1;printf: "T = %24.16e\n", T;printf: "n = %d\n", n;printf: "Nt = %d\n", Nt;printf: "N = %d\n", N;printf: "Data\n";printf: "t_k =";printf{k in 0..N}: " %24.16e\t", tau[k];printf: "\n";printf: "y =\n"; for {k in 0..N} {printf"k = %d\n", k;printf{j in 0..Nt, i in 1..n}: " %24.16e\n", y[k,j,i]; }end; # quit AMPL

In the following simulations, we take n=20 and N_t=400. By post-treatment of the results (see the *Scilab* file in the post IpOpt and AMPL use to solve time optimal control problems for an example), we obtain the following results:

- for $u^0=1$ and $u^1=5$: The minimal time obtained is $\mathtt{0.1689856}$.
- for $u^0=5$ and $u^1=1$: The minimal time obtained is $\mathtt{0.7267605}$.

In the above videos, red arrows correspond to Dirac masses in the control.

#### References

**[1]** J. Lohéac, E. Trélat, and E. Zuazua. *Control of the semi-discrete 1D heat equation under nonnegative control constraint.* In preparation.

**[2]** J. Lohéac, E. Trélat, and E. Zuazua. *Minimal controllability time for the heat equation under state constraints.* Mathematical Models and Methods in Applied Sciences, Volume 27, Issue 09, August 2017