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:

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