# Control of the semi-discrete 1D heat equation under nonnegative control constraint

##### 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 is the state to control the initial state and the Dirichlet control.
The continuous version of this problem has been analysed in [2] and it has been shown that for every initial state and every positive constant target , the controllability of this system under the control constraint

(1)

requires a positive waiting time as soon as the target state is different from the initial state .

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

(2)

where the matrices and are

(3)

where is the number of discretization points and (the component of ) stands for .
To obtain the above finite dimensional system, we have used centred finite differences.

The aim is then to minimize the time such that there exist a control steering the solution of (2) from to in time . For the sake of simplicity, we assume that and are constant vectors of , that is to say that

(4)

for some and Since the control shall satisfy the constraint (1), 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 .
Note that, using again the comparison principle, we can also show that if then, whatever the control is, the solution of (2) always satisfies .

To conclude, the optimization problem we aim to solve is

(5)

where and are given by (4) for some and , and where stands for the solution of (2) at time , with control and initial condition .

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 , the additional control constraint , and let increase to . 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 , defined by (5), there exists a non-negative control 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 Dirac masses.
In other words, the non-negative time optimal control takes the form:

with , and .
For a control of this form, the solution of (2) at time , with initial condition , is:

Consequently, the minimal time given by (5) is a minimizer of:

(6)

##### 3 Numerical implementation

The optimization problem (6) 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 , with and .
Let us then set for every .
We have and the optimization problem (6) also writes

In the above constraints, for every , where is the solution of (2), with control and initial condition . Notice that since is only constrained to be non-negative, and since the vector (given by (3)) is of the form , with , the constraint can be expressed as

Consequently, the parameters can be forgotten.

In order to numerically compute for , we are going to use the Crank-Nicholson method.
More precisely, given , we approximate by , with solution of

All in all, the fully discretized optimization problem, is

(7)

subject to the constraints

(8)

(9)

(10)

(11)

(12)

(13)

Note that instead of imposing , we chose to relax this condition in the constraint (12).
This choice has been a maid since we do not expect that the numerical solution exactly reaches the target . In practice, we will take, .

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

# define parameters
param n  = 20;             # number of spatial discretization points
param Nt = 400;            # number of time discretization points
param N  = floor((n+1)/2); # maximal number of control impulses
param eps= 1/(n*Nt);       # relaxation parameter  in (12)
param u0 = 1;              # parameters , see (4)
param u1 = 5;

# define variables
var y	  {k in 0..N, j in 0..Nt, i in 1..n};	# unknowns
var tau {k in 0..N} >=0;   # unknown impulse times  subject to (8)

# objective function, see (7)
minimize T:
sum {k in 0..N} tau[k];

# define the constraints
# constraint (13) for interior points
subject to y_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 (13) for i=n
subject to y_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 (13) for i=1
subject to y_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 (9)
subject to y_init {i in 1..n}:
y[0,0,i] = u0;
# terminal condition, constraint (12)
subject to y_end1 {i in 1..n}:
y[N,Nt,i] >= u1-eps;
subject to y_end2 {i in 1..n}:
y[N,Nt,i] <= u1+eps;
# constraint (10)
subject to continuity {k in 0..N-1, i in 1..n-1}:
y[k+1,0,i]=y[k,Nt,i];
# constraint (11)
subject to jump {k in 0..N-1}:
y[k+1,0,n]>=y[k,Nt,n];

# solve the problem with IpOpt
option solver ipopt;
option ipopt_options "max_iter=100000 linear_solver=mumps halt_on_ampl_error yes";
solve;

# display parameters and solution
printf: "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 and : The minimal time obtained is .

• for and : The minimal time obtained is .

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