Asymptotic method and transient terms in exact controls

Abstract

There is a narrow but hidden link between optimal control theory and the so-called Tikhonov regularization method. In fact, the small coefficient representing the marginal cost of the control can be interpreted as the regularization parameter in a Tikhonov method as far as there exists an exact control. This strategy enables one to adjust the cost function in the optimal control model in order to define the exact control which minimizes a given functional involving both the control but also the state variables during the control process. The goal of this paper is to suggest a method which gives a simple way to characterize and compute the exact control corresponding to the minimum of a given cost functional as said above. It appears as an extension of the phase control which is a finite dimensional version of the HUM control of J.L. Lions but for partial differential equations.

Keywords

Optimal control asymptotic methods Tikhonov regularization exact control dynamical systems

1. Introduction

Let us consider the classical finite dimensional model where A is a real square $N \times N$ matrix and B a rectangular one with dimensions $N \times m$ ( $m ⩽ N$ ). Both could be time dependent. The state variable is denoted by X and it is a vector of $R^{N}$ depending on the time t, solution of: $\begin{matrix} (1) & \frac{d X}{d t} = A X + F + B u, X (0) = X_{0}, \end{matrix}$ where $u = {u_{i}} \in L^{2} (] 0, T [; R^{m}) = [L^{2} (] 0, T {[)]}^{m}$ is the control, $F \in [L^{2} (] 0, T {[)]}^{N}$ , F is a given right-hand side, and $X_{0} \in R^{N}$ is the initial condition. The existence and uniqueness of $X \in [H^{1} (] 0, T {[)]}^{N}$ are standard. We introduce a first guess of a path from $X_{0}$ to $X_{d} (T)$ – say $X_{d} (t)$ – such that $X_{d} \in [H^{1} (] 0, T {[)]}^{N}$ . For instance, in a routing process, $X_{d}$ could be a solution obtained by a coarse algorithm or a first guess of a path which seems to be the most interesting by the user but not necessarily neither optimal nor feasible. For instance, for a sailing boat, the path face to the wind can be the shortest path in terms of distance but impossible to follow.

Then we define a criterion by ( $(\cdot, \cdot)$ is the scalar product in $R^{N}$ or $R^{m}$ , and the norm of a vector – say v – is $| v |$ ; finally the composition between matrix and vectors is represented by a dot: . furthermore the term $d t$ will be omitted when there is no ambiguity): $\begin{matrix} (2) & \begin{matrix} J^{ε} (u) & = \frac{1}{2} (D . (X (T) - X d (T)), X (T) - X d (T)) \\ + \frac{ε}{2} \int_{0}^{T} {(C . (X (t) - X_{d} (t)), X (t) - X_{d} (t)) + (R . u (t), u (t))} (t) d t, \end{matrix} \end{matrix}$ where D, C (respectively R) are $N \times N$ matrices (respectively a $m \times m$ matrix). They are real symmetrical and positive matrices characterizing the optimal control model that we consider in this paper. The two matrices C and D are not necessarily definite but R is assumed to be always definite. Their definitions are the decision of the engineers who plan to use the optimal control algorithm. One can also choose for C a piecewise continuous function of time and for instance it can be zero outside of an interval $] t_{1}, t_{2} [$ . In fact C (but also D) can imply only few components of $R^{N}$ . The matrix D is definite as far as all the components of $X_{d} (T)$ should be reached at time T. This the case in the following in order to simplify the presentation. Finally, the optimal control model considered, consists in solving the following optimization problem which is in fact a so-called optimal control problem: $\begin{matrix} (3) & min_{v \in L^{2} (] 0, T [; R^{m})} J^{ε} (v) . \end{matrix}$ One could also discussed the possibility to take into account a limitation on the control variable. For instance one could consider the following admissible set for the control (where $c_{i} > 0$ ): $\begin{matrix} (4) & U_{ad} = {v = {v_{i}} \in L^{2} (] 0, T [; R^{m}), \forall i = 1, m : | v_{i} | ⩽ c_{i}} . \end{matrix}$ The basic point is to use a non empty closed and convex set of the space $L^{2} (] 0, T [; R^{m})$ . As far as one wishes to handle a more smooth control ( $C^{0} ([0, T]; R^{m})$ for instance), one can introduce the norm $H^{1} (] 0, T [; R^{m})$ for the control in the definition of the criterion instead of the one in $L^{2} (] 0, T [; R^{m})$ . If so, one can introduce an alternative criterion to be minimized by setting (for instance): $\begin{matrix} (5) & H^{ε} (v) = J^{ε} (v) + \frac{ε}{2} \int_{0}^{T} (R_{d} . \frac{d u}{d t}, \frac{d u}{d t}) (t) d t, \end{matrix}$ where $R_{d}$ is another symmetrical and positive matrix. In fact it should be definite for the component of the control u that are required to be $C^{0} ([0, T])$ . Furthermore one can prescribe boundary conditions (at $t = 0$ and $t = T$ ) for the components of u which are in the space $H^{1} (] 0, T [)$ . The extension ot the present analysis given in this paper to this criterion (5) is straight forward excepted may be in the expression of the control with respect to the adjoint state that we made explicit at Remarque 1. Hence we do not focus too much on it in the following for sake of simplicity mainly in the notations.

The existence and uniqueness of a solution to the optimization model (3) or (5), are classical because of the linearity of the state equation and the strict convexity of the criterion with respect to the control variable u. Let us denote it by $u^{ε}$ . Our goal is to study the asymptotic behavior of this solution when $ε \to 0$ . It is proved in this paper, that the limit exists iff the model (1) can be exactly controlled at time T. In other words, it is necessary that there exists a control $v \in [L^{2} (] 0, T {[)]}^{m}$ such that the solution of (1) obtained with the control v, satisfies: $X (T) = X_{d} (T)$ . Furthermore, this limit is the unique one which minimizes the term in factor of ε among the set of all the exact controls. This is the important result of the paper. Because this property is very close to the Tikhonov method for solving linear systems associated to a symmetrical and positive operator but not definite we mention that it is a Tikhonov regularization.

2. Characterization of $u^{ε}$

The optimality conditions of problem (3) are the following ones ( $^{tr}$ denotes the transposed of a matrix): $\begin{matrix} (6) & \{\begin{matrix} \forall t \in [0, T], \frac{d X^{ε}}{d t} = A . X^{ε} + F + B . u^{ε}, X^{ε} (0) = X_{0}, \\ \forall t \in [0, T], \frac{d P^{ε}}{d t} = -^{tr} A . P^{ε} - ε C . (X^{ε} - X_{d}), P^{ε} (T) = D . (X^{ε} (T) - X_{d} (T)), \\ \forall t \in [0, T],^{tr} B . P^{ε} + ε R . u^{ε} = 0 . \end{matrix} \end{matrix}$ One can refer to R. Bellman [2,3], L. Pontryaguin [19], J.L. Lions [14] or P. Faurre and M. Robin [8]. Concerning a discretized formulation in time the classical references are those of R. Kalman [11–13].

Remark 1.
In the case where the criterion $H^{ε}$ (see (5)) would be used, the optimality equation (the third equation of system (6) becomes (we do not prescribe the values of the control at $t = 0$ and $t = T$ and therefore we get homogeneous Neumann boundary conditions): $\begin{matrix} (7) & ^{tr} B . P^{ε} + ε (R . u^{ε} - R_{d} . \frac{d^{2} u^{ε}}{d t^{2}}) = 0, and \frac{d u^{ε}}{d t} (0) = \frac{d u^{ε}}{d t} (0) = 0 . \end{matrix}$ If the boundary conditions on $u^{ε}$ are prescribed at $t = 0$ or/and $t = T$ they would replace those explicited at (7).
Remark 2.
There are at least three methods for solving the optimality equations (6). One which is not realistic because of the volume of the computation for transient problem, consists in solving globally the linear system (6). Another possibility, which is the most used in modern application, is to use a gradient algorithm with a preconditioning or even a conjugate gradient also with a preconditioning. A basic advantage is to be able to treat easily the constraints on the control if there are some. But the robustness of the solution method always depends on the fact that a controllability property of the system is satisfied. When the parameter ε is small (see [6]) this controllability property is necessary. If it is not satisfied the solution exploses when $ε \to 0$ . A third strategy which was very popular among automaticians is to use the Riccati method [2–12]. Unfortunately it can’t always be applied. Let us consider a case where it can be used. We assume that we know a first guess of a control law – say $u_{d}$ – and $X_{d}$ which is solution of the state equation: $\begin{matrix} (8) & \frac{d X_{d}}{d t} = A . X_{d} + F + B . u_{d}, X_{d} (0) arbitrary . \end{matrix}$ In fact if we prescribe (recommended) $X_{d} (0) = X_{0}$ , $X_{d}$ is the solution of the state equation with a control equal to $u_{d}$ . Let us set: $\begin{matrix} (9) & u = u_{d} + δ u, \end{matrix}$ and the new control is $δ u$ . In this case $X_{d}$ should satisfies: $\begin{matrix} (10) & \frac{d X_{d}}{d t} = A . X_{d} + F + B . u_{d}, X_{d} (0) = X_{0} . \end{matrix}$ Because the new control is $δ u$ it is the one which appear in the expression of the cost functional. Therefore the optimality equation can be written: $\begin{matrix} ^{tr} B . P^{ε} + ε R . δ u^{ε} = 0 . \end{matrix}$ The first step in Riccati method, consists in making a change of variables setting: $\begin{matrix} (11) & Y (t) = X^{ε} (t) - X_{d} (t), Y (0) = X^{ε} (0) - X_{d} (0) = X_{0} - X_{d} (0) . \end{matrix}$ This strategy consists in looking for a matrix – say Z – which is symmetrical, time dependent and such that (the control involved in the expression of the cost functional is $δ u^{ε}$ ): $\begin{matrix} \forall t \in [0, T], P^{ε} (t) = Z (t) . (X^{ε} (t) - X_{d} (t)) = Z (t) . Y (t) and^{tr} B . P^{ε} = - ε R . δ u^{ε} . \end{matrix}$ By substituting this assumed expression of $P^{ε}$ into the equations (6), one obtains a necessary relation which should be satisfied by Z: $\begin{matrix} (12) & \begin{matrix} \frac{d Z}{d t} . Y + Z . \frac{d Y}{d t} = -^{tr} A . Z . Y - ε C . Y and \frac{d Y}{d t} = A . Y - \frac{1}{ε} R^{- 1} .^{tr} B . Z . Y, \\ with \\ Z (T) . Y (T) = D . Y (T), Y (0) = X_{0} - X_{d} (0) . \end{matrix} \end{matrix}$ Or else: $\begin{matrix} (13) & [\frac{d Z}{d t} + Z . A +^{tr} A . Z + ε C - \frac{1}{ε} Z . B . R^{- 1} .^{tr} B . Z] . Y = 0 and Z (T) . Y (T) = D . Y (T), \end{matrix}$ and finally, if Z is solution of the so-called Riccati equation: $\begin{matrix} (14) & \frac{d Z}{d t} + Z . A +^{tr} A . Z + ε C - \frac{1}{ε} B . R^{- 1} .^{tr} B . Z = 0 and Z (T) = D \end{matrix}$ Finally we have found the solution of the optimality equation. Practically the Riccati equation can be solved using a Runge–Kutta scheme but it is necessary for the stability to use at least a one step of a Newton algorithm for the non linear terms. If there are bounds on the control, the method is no more rigorous. In this case, people very often use an approximation of the optimal control by setting ( $P_{U_{ad}}$ is the euclidian projection on the admissible set for the controls – say $U_{ad}$ – assumed to be convex, non empty and containing 0 as an internal point): $\begin{matrix} (15) & u^{ε} (t) = P_{U_{ad}} [- \frac{1}{ε} R^{- 1} .^{tr} B . Z . (X^{ε} - X_{d})] . \end{matrix}$

In fact we suggest hereafter, a different method for solving the optimality equation based on the phase control strategy (terminology used in 1D by automation engineers) which is obtained as the limit of the optimal control when $ε \to 0$ assuming the existence of an exact control. In partial differential equations this method has been introduced differently by J.L. Lions [16] and is the so-called HUM method. But usually the authors of papers concerning exact control papers didn’t define the phase control as the limit of the optimal control (when the exact controllability assumption is satisfied) when the cost of the control tends to zero. And our goal in the following is precisely to study the asymptotic behavior of $u^{ε}$ when $ε \to 0$ . For this purpose we make use of the asymptotic method following the method used in [6] for optimal control model.
2.1. The formal asymptotic expansion

Our goal in this section is to construct an asymptotic expansion with respect to ε of the optimal control and of the state function associated. The first step is formal (this subsection and the mathematical justification is presented in the second step; see Section 2.2). Let us set a priori (see J.L. Lions [15]): $\begin{matrix} (16) & \{\begin{matrix} X^{ε} = X^{0} + ε X^{1} + \dots \\ P^{ε} = P^{0} + ε P^{1} + \dots \\ u^{ε} = u^{0} + ε u^{1} + \dots \end{matrix} \end{matrix}$ By introducing these expressions in the equations (6) and by equating the terms with the same power in ε, one gets (let us recall that C can depend on the time variable): $\begin{matrix} (17) & \begin{matrix} Order 0 \{\begin{matrix} \forall t \in [0, T], \frac{d X^{0}}{d t} = A . X^{0} + B . u^{0} + F, X^{0} (0) = X_{0}, \\ \forall t \in [0, T], \frac{d P^{0}}{d t} = -^{tr} A . P^{0}, P^{0} (T) = D . (X^{0} (T) - X_{d} (T)), \\ \forall t \in [0, T],^{tr} B . P^{0} = 0, \end{matrix} \\ Order 1 \{\begin{matrix} \frac{d X^{1}}{d t} = A . X^{1} + B . u^{1}, X^{1} (0) = 0, \\ \frac{d P^{1}}{d t} = -^{tr} A . P^{1} - C . (X^{0} - X_{d}), P^{1} (T) = D . X^{1} (T), \\ ^{tr} B . P^{1} + R . u^{0} = 0, \end{matrix} \\ and so on \dots \end{matrix} \end{matrix}$ Let us sketch a solution method for computing a few of the terms of this assumed asymptotic expansion. First of all, we require a standard controllability assumption given in the next statement.

Assumption 1 (Controllability (see R. Bellman [2] or I. Pontryagin [19])).

Let us consider a vector $P \in {[C^{0} ([0, T])]}^{N})$ depending on time t and satisfying the relations: $\begin{matrix} \forall t \in [0, T] :^{tr} B . P = 0 and \frac{d P}{d t} = -^{tr} A . P \end{matrix}$ Then the system ( 1 ) is said to be exactly controllable iff $P = 0$ .

(The terminology is explained in the following; see ( 18 )).

Assuming the exact controllability property recalled in Assumption 1, one can claim from (17) that a necessary condition for this assumed asymptotic expansion, one should have $P^{0} = 0$ . This also implies that: $\begin{matrix} (18) & D . (X^{0} (T) - X_{d} (T)) = 0, \end{matrix}$ and if D is definite (which is the case considered in this text) one deduces that the state $X_{d} (T)$ can be reached at time T; it is the exact controllability property as soon as one is able to compute the corresponding control $u^{0}$ . This is the goal that we try to achieve hereafter by a similar method to the one introduced in [6].

Remark 3.
In fact the controllability Assumption as formulated previously, is certainly too restrictive if the matrix D is singular. But it is the one used in many books on the subject. Furthermore, if one wishes to introduce a sufficient condition for the controllability of $D . X^{0} (T)$ which is weaker than the one of $X^{0} (T)$ , the compatibility between the three matrices A, B and D would be a little bit more difficult to handle technically, even if this would be possible.

The control $u^{0}$ should be given from $P^{1}$ by the following expression (this is a necessary condition but not a sufficient one): $\begin{matrix} (19) & u^{0} = - R^{- 1} .^{tr} B . P^{1}, \end{matrix}$ and $P^{1}$ should be solution of: $\begin{matrix} (20) & \begin{matrix} \frac{d P^{1}}{d t} = -^{tr} A . P^{1} - C . (X^{0} - X_{d}) with P^{1} (0) = Φ \in R^{N} \\ where Φ and X^{0} are still unknown at this step . \end{matrix} \end{matrix}$ Let us point out that $X^{0}$ is unknown because it depends on the control $u^{0}$ which is still undetermined. One has a closed loop problem as far as $u^{0}$ and $X^{0}$ are both unknowns at this step and it is not obvious to eliminate one of these two with respect to the other one. It is the goal of the next subsection to solve this difficulty due to the transient term as far as $C \neq 0$ . And this is the main point treated in this paper.
2.2. Computation of $u^{0}$ and $X^{0}$

There is a constraint on $X^{0}$ due to equation (18). But it also depends on $u^{0}$ which depends on $P^{1}$ solution (20) where the term $X^{0}$ also appears through the symmetrical and positive (but not necessarily definite) matrix C which furthermore can depend on the time variable. It can be zero for some time intervals and even can have a kernel where it is not identically zero. At this step, it will only be possible to compute the components of $X^{0}$ which are orthogonal to the kernel of C (in the range of C because C is always symmetrical) and therefore for any $t \in [0, T]$ : $\begin{matrix} R^{N} = Ker (^{tr} C (t)) \oplus Range (C (t)) = Ker (C (t)) \oplus Range (C (t)) . \end{matrix}$ In order to avoid any misinterpretation we denote by $X_{C}^{0}$ this component of $X^{0}$ which is in the range of C and thus one should find $X^{0}$ such that $C . X^{0} = C . X_{C}^{0}$ . Obviously if $C (t)$ is definite any possible confusion disappears. For sake of clarity, we introduce the functional space in which we look for $X_{C}^{0}$ by (a.e. t means: almost every t): $\begin{matrix} (21) & V_{C} = {a.e. t \in] 0, T [; δ C (t) \in Range (C (t)) and δ C \in [L^{2} (] 0, T [)]^{N}} . \end{matrix}$ For any $δ X^{0} \in V_{C}$ let us now introduce a vector field Q solution of: $\begin{matrix} (22) & \frac{d Q}{d t} = -^{tr} A . Q - C . δ X^{0} with Q (0) = δ Φ \in R^{N} . \end{matrix}$ By multiplying the first equation in (17) by Q and after several integrations by parts, one obtains (let us point out that from the definition of $X_{C}^{0}$ one has: $C . X^{0} = C . X_{C}^{0}$ ): $\begin{matrix} (23) & \{\begin{matrix} \forall (δ Φ, δ X^{0}) \in R^{N} \times V_{C}, \\ \int_{0}^{T} (R^{- 1} .^{tr} B . P^{1} (t),^{tr} B . Q (t)) d t + \int_{0}^{T} (C . X_{C}^{0}, δ X^{0}) \\ = \int_{0}^{T} (F, Q) + (X_{0}, δ Φ) - (X (T), Q (T)) . \end{matrix} \end{matrix}$ Assuming for instance that D is definite, one has $X (T) = X_{d} (T)$ and therefore: $\begin{matrix} (24) & \{\begin{matrix} \forall (δ Φ, δ X^{0}) \in R^{N} \times V_{C}, \\ \int_{0}^{T} ({R^{- 1}}^{tr} B . P^{1} (t),^{tr} B . Q (t)) d t + \int_{0}^{T} (C . X_{C}^{0}, δ X^{0}) \\ = \int_{0}^{T} (F, Q) + (X_{0}, δ Φ) - (X_{d} (T), Q (T)) . \end{matrix} \end{matrix}$ It is worth noting at this step that $P^{1}$ , $X_{C}^{0}$ depends on both Φ, $X_{0}$ and also on $X_{d}$ which appears as a data in equation (24). In order to solve (24) we split $P^{1}$ into the sum of three terms:

$P_{d}^{1}$ solution of the following equation and which only depends on $X_{d}$ which is known and therefore can be computed separately: $\begin{matrix} (25) & \frac{d P_{d}^{1}}{d t} = -^{tr} A . P_{d}^{1} + C . X_{d}, P_{d}^{1} (0) = 0 . \end{matrix}$

$P_{x^{0}}^{1}$ (to be determine as a function of $X_{C}^{0}$ which is unknown at this step) solution of: $\begin{matrix} (26) & \frac{d P_{x^{0}}^{1}}{d t} = -^{tr} A . P_{x^{0}}^{1} - C . X_{C}^{0}, P_{x^{0}}^{1} (0) = 0 . \end{matrix}$

${\tilde{P}}^{1}$ (dependent only on Φ to be determined), solution of: $\begin{matrix} (27) & \frac{d {\tilde{P}}^{1}}{d t} = -^{tr} A . {\tilde{P}}^{1}, {\tilde{P}}^{1} (0) = Φ . \end{matrix}$

One has to determine both Φ and

X_{C}^{0}

R^{N} \times V_{C}

such that (let us point out that at this step we only look for

X_{C}^{0}

in the space

V_{C}

and Q being solution of (22) and we assume for sake of simplicity that D which appears in the criterion to be minimized as an observer operator of the final condition, is definite):

\begin{matrix} (28) & \{\begin{matrix} \forall (δ Φ, δ X^{0}) \in R^{N} \times V_{C}, \\ \int_{0}^{T} (R^{- 1} .^{tr} B . {\tilde{P}}^{1},^{tr} B . Q) + \int_{0}^{T} (R^{- 1} .^{tr} B P_{x^{0}}^{1},^{tr} B . Q) + \int_{0}^{T} (C . X^{0}, δ X^{0}) \\ = \int_{0}^{T} (F, Q) + (X_{0}, δ Φ) - (X_{d} (T), Q (T)) - \int_{0}^{T} (R^{- 1} .^{tr} B . P_{d}^{1},^{tr} B . Q) \end{matrix} \end{matrix}

We shall check that

X_{C}^{0}

is a part of the solution of (1) if

u = u^{0}

given by (19). This is not straightforwards for at least three reasons:

first because $X_{C}^{0}$ is looked for in the space $[L^{2} (] 0, T {[)]}^{N}$ , therefore a regularity result has to be proved;

secondly, because the initial ( $t = 0$ ) and final ( $t = T$ ) conditions should be satisfied by $X_{C}^{0}$ ;

thirdly, because $X_{C}^{0}$ should satisfy (at least partially) the equation (1) with the control $u^{0}$ defined by (19), once $P^{1}$ will be defined.

We are going to use Lax-Milgram Theorem (see J. Cea for instance [4] for solving the equation (23) or equivalently (28). Let us first introduce few notations:

\begin{matrix} (29) & \{\begin{matrix} Z = (Φ, X_{C}^{0}) \in R^{N} \times V_{C} δ Z = (δ Φ, δ X^{0}) \in R^{N} \times V_{C}, \\ Λ (Z, δ Z) = (C . X_{C}^{0}, δ X^{0}) + \int_{0}^{T} [({R^{- 1}}^{tr} B . {\tilde{P}}^{1},^{tr} B . Q) + (^{tr} B . P_{x^{0}}^{1},^{tr} B . Q)] \\ where {\tilde{P}}^{1} and Q should respectively satisfies (27) and (22): \\ L (δ Z) = (X_{0}, δ Φ) - (X_{d} (T), Q (T)) + \int_{0}^{T} [(F, Q) - ({R^{- 1}}^{tr} B P_{d}^{1},^{tr} B . Q)] . \end{matrix} \end{matrix}

The problem to be solved is therefore: $\begin{matrix} (30) & \{\begin{matrix} Find Z = (Φ, X_{C}^{0}) \in R^{N} \times V_{C} such that: \\ \forall δ Z \in R^{N} \times V_{C}, Λ (Z, δ Z) = L (δ Z) . \end{matrix} \end{matrix}$ One has the following properties:

L is a linear and continuous form on the Hilbert space: $R^{N} \times V_{C}$ ;

Λ is a bilinear positive and symmetrical form on: $R^{N} \times V_{C}$ ;

if Z satisfies $Λ (Z, Z) = 0$ one has $X_{C}^{0} = 0$ and $^{tr} B ({\tilde{P}}^{1} + P_{x^{0}}^{1}) = 0$ . But if $X_{C}^{0} = 0$ then $P_{x^{0}}^{1} = 0$ and therefore $^{tr} B {\tilde{P}}^{1} = 0$ . Because of (26) and from the controllability assumption (see Assumption 1), one gets ${\tilde{P}}^{1} = 0$ . Finally, Λ is positive definite on $R^{N} \times V_{C}$ .

The last point consists in proving the coerciveness of Λ on the space

R^{N} \times V_{C}

(the difficulty is due to the fact the

V_{C}

is not necessarily a finite dimensional space because its elements depend on the time variable).

Let us notice that for any $0 < α < 1$ (using the triangular Cauchy–Schwarz inequality also called Yung inequality1

¹
$\forall a, b \in R$ , $α \in R^{+} : 2 a b ⩽ α a^{2} + \frac{1}{α} b^{2}$ .

) and denoting by

r_{0}

(respectively by

r_{1}

) the smallest (respectively largest) eigenvalue of the

m \times m

matrix R (and their inverses for

R^{- 1}

with obviously

r_{0} ⩽ r_{1}

\begin{matrix} (31) & \{\begin{matrix} Λ (Z, Z) = \int_{0}^{T} (R^{- 1} [^{tr} B . ({\tilde{P}}^{1} + P_{x^{0}}^{1})],^{tr} B . ({\tilde{P}}^{1} + P_{x^{0}}^{1})) + (C . X_{C}^{0}, X_{C}^{0}) \\ ⩾ \int_{0}^{T} (C . X_{C}^{0}, X_{C}^{0}) + \int_{0}^{T} [\frac{1}{r_{1}} (1 - α) |^{tr} B {\tilde{P}}^{1} |^{2} + \frac{1}{r_{0}} (1 - \frac{1}{α}) |^{tr} B . P_{x^{0}}^{1} |^{2}] . \end{matrix} \end{matrix}

But from the equation (26) which defines

P_{x^{0}}^{1}

, one gets:

\begin{matrix} (32) & P_{x^{0}}^{1} (t) = - \int_{0}^{t} e^{^{tr} A (s - t)} . C . X_{C}^{0} (s) d s, \end{matrix}

which implies that there exists a constant

c_{0}

depending on T and on the matrices A, B and C such that:

\begin{matrix} (33) & \int_{0}^{T} {|^{t} B . P_{x^{0}}^{1} (s) |}^{2} d s ⩽ c_{0} \int_{0}^{T} (C . X_{C}^{0} (s), X_{C}^{0} (s)) d s . \end{matrix}

Because

1 - \frac{1}{α} < 0

and choosing α such that:

\frac{c_{0}}{c_{0} + r_{0}} < α < 1

(which is always possible), one obtains that there exists a strictly positive constant

c_{1}

independent on Z and noting that C is positively definite on

V_{C}

, one gets (

C . X_{C}^{0} = X_{C}^{0}

\begin{matrix} (34) & \{\begin{matrix} \forall Z = (Φ, X^{0}) \in R^{N} \times V_{C}, \\ Λ (Z, Z) ⩾ c_{1} [\int_{0}^{T} | X_{C}^{0} (s) |^{2} d s + \int_{0}^{T} |^{tr} B . {\tilde{P}}^{1} (s) |^{2} d s], \end{matrix} \end{matrix}

and finally because of the exact controllability Assumption 1 which ensures that:

\begin{matrix} (35) & Φ \in R^{N} \to \sqrt{\int_{0}^{T} {|^{tr} B . {\tilde{P}}^{1} (s) |}^{2} d s}, \end{matrix}

is a norm on

R^{N}

(equivalent to any other one because we are in a finite dimensional space), there exists a constant

c_{2} > 0

independent on both Φ and

X^{0}

such that:

\begin{matrix} (36) & Λ (Z, Z) ⩾ c_{2} [\int_{0}^{T} {| X_{C}^{0} (s) |}^{2} d s + | Φ |^{2}], \end{matrix}

which establishes the coerciveness of Λ on the space

R^{N} {\times V}_{C}

. Therefore, the Lax-Milgram Theorem enables one to claim that there is a unique solution to the equation (30).

Let us summarized and complete the results obtained at this step in the following statement.

Theorem 2.1.

Let us denote by $(Z^{0} = (Φ^{0}), X_{C}^{0} \in R^{N} \times V_{C}$ the unique solution of ( 30 ). We set: $\begin{matrix} (37) & u^{0} = -^{tr} B P^{1}, \end{matrix}$ where $P^{1}$ is the solution of ( 24 ). Then $u^{0}$ is an exact control for the initial condition $X_{0} \in R^{N}$ and the guessed path $X_{d}$ (i.e. $X^{0} (T) = X_{d} (T)$ ).

Proof.

Let us introduce Q solution of: $\begin{matrix} (38) & \dot{Q} = -^{tr} A . Q, Q (0) = δ Φ \in R^{N} . \end{matrix}$ By multiplying by Q the state equation (1) satisfied by $X^{0}$ associated to the control $u^{0}$ , and because of the relation satisfied by $P^{1}$ , one obtains: $\begin{matrix} (39) & \forall δ Φ \in R^{N} (X^{0} (T), Q (T)) = (X_{d} (T), Q (T)) . \end{matrix}$ In fact by choosing $\begin{matrix} δ Φ = (e^{^{tr} A T}) . (X^{0} (T) - X_{d} (T)) \end{matrix}$ which implies that $\begin{matrix} Q (T) = X^{0} (T) - X_{d} (T) \end{matrix}$ one can conclude that: $\begin{matrix} X^{0} (T) = X_{d} (T) . \end{matrix}$ □

Furthermore, the solution $X_{C}^{0}$ is more regular than the framework of the Lax-Milgram Theorem suggests. This is given in the following result.

Theorem 2.2.

The component $X_{C}^{0}$ of the solution of equation ( 30 ) which is in the space $[L^{2} (] 0, T {[)]}^{N}$ , is more regular and is at least in the space $[H^{1} (] 0, T {[)]}^{N}$ .

Proof.

From the equation (30) one has (choosing $δ Φ = 0$ in (30)): $\begin{matrix} (40) & \{\begin{matrix} \forall δ X^{0} \in V_{C} : \\ \int_{0}^{T} (C . X_{C}^{0}, δ X^{0}) + \int_{0}^{T} (^{tr} B . P_{x^{0}}^{1},^{tr} B . Q) + \int_{0}^{T} (^{tr} B . P_{d}^{1},^{tr} B . Q) \\ + \int_{0}^{T} (^{tr} B . {\tilde{P}}^{1},^{tr} B . Q) = - (X_{d}, Q (T)), \\ where Q is solution of (let us recall that C can be a function of time): \\ \frac{d Q}{d t} = -^{tr} A . Q - C . δ X^{0}, Q (0) = 0 . \end{matrix} \end{matrix}$ The solution for Q is: $\begin{matrix} (41) & Q (t) = - \int_{0}^{t} e^{^{tr} A (s - t)} . C (s) . δ X^{0} (s) d s \end{matrix}$ Therefore for any function $G \in [L^{2} (] 0, T {[)]}^{N}$ one has from Fubini Theorem [5]: $\begin{matrix} (42) & \{\begin{matrix} \int_{0}^{T} (^{tr} B . G,^{tr} B . Q) = - \int_{0}^{T} (e^{^{tr} A s} . C (s) . δ X^{0} (s), [\int_{s}^{T} e^{- A t} . B .^{tr} B . G (t) d t]) d s \\ = - \int_{0}^{T} (δ X^{0} (s), \int_{s}^{T} [C (s) . e^{A (s - t)} . B .^{tr} B . G (t) d t]) d s . \end{matrix} \end{matrix}$ Finally we proved that (40) leads to (let us recall that C is a symmetrical and positive matrix): $\begin{matrix} (43) & \{\begin{matrix} C (s) . X_{C}^{0} (s) = C (s) . e^{A (s - T)} . X_{d} - C (s) \int_{s}^{T} e^{A (s - t)} . B .^{tr} B . ({\tilde{P}}^{1} + P_{d}^{1} + P_{x^{0}}^{1}) (t) d t, \\ but C has an inverse on its range, say C_{R a}^{- 1} and X_{c}^{0} belongs to it: \\ X_{C}^{0} (s) = C_{R a}^{- 1} (s) . C . [e^{A (s - T)} . X_{d} - \int_{s}^{T} e^{A (s - t)} . B .^{tr} B . ({\tilde{P}}^{1} + P_{d}^{1} + P_{x^{0}}^{1}) (t) d t] . \end{matrix} \end{matrix}$ Because ${\tilde{P}}^{1}$ , $P_{d}^{1}$ and $P_{x^{0}}^{1}$ are at least $[H^{1} (] 0, T {[)]}^{N}$ one can claim that $X^{0}$ has also the same regularity where C is definite and therefore admits an inverse. If C has an inverse on all the segment $[0, T]$ then $X_{C}^{0}$ belongs to the space $[H^{1} (] 0, T {[)]}^{N}$ and in this particular case one has: $\begin{matrix} (44) & X_{C}^{0} (s) = X^{0} = e^{A (s - T)} . X_{d} - \int_{s}^{T} e^{A (s - t)} . B .^{tr} B . ({\tilde{P}}^{1} + P_{d}^{1} + P_{x^{0}}^{1}) (t) d t . \end{matrix}$ □

Conversely, the term $X_{C}^{0}$ characterized at (43) is in fact locally the component of $X^{0}$ in the range of C solution of (1) with the control $u^{0}$ . If C is invertible everywhere in $[0, T]$ then $X_{C}^{0} = X^{0}$ and satisfies the final condition $X_{C}^{0} (T) = X^{0} (T)$ . These results are summarized in the following statement.

Theorem 2.3.

Let $X^{0}$ the solution of ( 1 ) with the control: $\begin{matrix} (45) & u^{0} = -^{tr} B . P^{1} = -^{tr} B . ({\tilde{P}}^{1} + P_{x^{0}}^{1} + P_{d}^{1}) . \end{matrix}$ From Theorem 2.1 one has $X^{0} (T) = X_{d}$ . But in addition, if $X_{C}^{0}$ is solution of ( 23 ), one has also: $\begin{matrix} (46) & C . (X_{C}^{0} - X^{0}) = 0 . \end{matrix}$ If C is invertible everywhere over $[0, T]$ then: $\begin{matrix} (47) & X^{0} = X_{C}^{0} everywhere over] 0, T [. \end{matrix}$ For instance, if C is invertible in $] t_{1}, t_{2} [$ then: $\begin{matrix} (48) & X^{0} = X_{C}^{0} but only over this time interval . \end{matrix}$ If $C = 0$ in $] t_{1}, t_{2} [$ then $X_{C}^{0} = 0$ on $] t_{1}, t_{2} [$ .

Proof.

Let us start from the equation (1) and let us introduce in it, the expression of $u^{0}$ given in the formulation of Theorem 2.3. This leads to (Q is solution of (22)): $\begin{matrix} (49) & \{\begin{matrix} \forall (δ Φ, δ X^{0}) \in R^{N} \times [L^{2} (] 0, T {[)]}^{N} : \\ \int_{0}^{T} (\frac{d X^{0}}{d t} - A . X^{0} - F + B .^{tr} B . P^{1}), Q) = 0, X^{0} (0) = X_{0} . \end{matrix} \end{matrix}$ From an integration by parts, one gets: $\begin{matrix} (50) & (X^{0} (T), Q (T)) - (X_{0}, δ Φ) - \int_{0}^{T} (X^{0}, \frac{d Q}{d t} +^{tr} A . Q) + \int_{0}^{T} (^{tr} B . P^{1},^{tr} B . Q) = 0, \end{matrix}$ and from the definition of Q given at (22) (pay attention that it is $X^{0}$ and not $X_{C}^{0}$ which appears in the following equation): $\begin{matrix} (51) & (X^{0} (T), Q (T)) - (X_{0}, δ Φ) + \int_{0}^{T} (C . X^{0}, δ X^{0}) + \int_{0}^{T} (^{tr} B . P^{1},^{tr} B . Q) = 0 . \end{matrix}$ Finally, using the definition of $P^{1}$ at (24): $\begin{matrix} (52) & \forall (δ Φ, δ X^{0}) \in R^{N} \times V_{C} : (X^{0} (T) - X_{d} (T), Q (T)) + \int_{0}^{T} (C . (X^{0} - X_{C}^{0}), δ X^{0}) = 0 . \end{matrix}$ One can choose $δ X_{C}^{0}$ arbitrary in (52). Because of Theorem 2.1, this implies that: $\begin{matrix} C . (X^{0} - X_{C}^{0}) = 0 . \end{matrix}$ Thus Theorem 2.3 is proved. □

The next step consists in proving a convergence result for $u^{ε}$ to $u^{0}$ when $ε \to 0$ and the optimality for $u^{0}$ with respect to the norm chosen in the cost criterion for the control (Tikhonov property [20,21]).

2.3. Convergence result when

ε \to 0

Let us choose $u = u^{0}$ in the criterion $J^{ε}$ defined at (2). This leads to (because $u^{0}$ is exact i.e. $X^{0} (T) = X_{d}$ as far as D is still assumed to be definite): $\begin{matrix} (53) & \{\begin{matrix} \frac{ε}{2} \int_{0}^{T} [| u^{ε} |^{2} (t) + (C . (X^{ε} - X_{d}), (X^{ε} - X_{d})) (t)] d t \\ + \frac{1}{2} (D . (X^{ε} (T) - X_{d} (T)), (X^{ε} (T) - X_{d} (T))) \\ = J^{ε} (u^{ε}) ⩽ J^{ε} (u^{0}) = \frac{ε}{2} \int_{0}^{T} [| u^{0} |^{2} (t) + (C . (X^{0} - X_{d}), (X^{0} - X_{d})) (t)] d t . \end{matrix} \end{matrix}$ This implies that the sequence $u^{ε}$ is bounded in the space $[L^{2} {(] 0, T [)}^{N}$ . From the equation satisfied by $X^{ε}$ one can claim that $X^{ε}$ is also bounded in the space $[H^{1} (] 0, T {[)]}^{N}$ and from (53) that: $\begin{matrix} (54) & lim_{ε \to 0} X^{ε} (T) = X_{d} (T) . \end{matrix}$ From the previous statement, one can extract from $(u^{ε}, X^{ε})$ a subsequence denoted by $(u^{ε^{'}}, X^{ε^{'}})$ and such that: $\begin{matrix} (55) & \{\begin{matrix} u^{ε^{'}} ⇀ u^{*} weakly in [L^{2} (] 0, T {[)]}^{m}, \\ X^{ε^{'}} ⇀ X^{*} weakly in [H^{1} (] 0, T {[)]}^{N} \subset C^{0} {([0, T])}^{N} . \end{matrix} \end{matrix}$ Furthermore, one has from the inequality (53): $\begin{matrix} (56) & X^{*} (T) = X_{d} (T) in R^{N} . \end{matrix}$ From the equations satisfied by the sequence $(u^{ε}, X^{ε})$ and the semi-lower continuity of the $L^{2}$ -norm with respect to the weak topology, one obtains by taking the limit: $\begin{matrix} (57) & \{\begin{matrix} \frac{d X^{*}}{d t} = A . X^{*} + B . u^{*}, X^{*} (0) = X_{0}, \\ \int_{0}^{T} [| u^{*} |^{2} (t) + (C . (X^{*} - X_{d}), (X^{*} - X_{d})) (t)] d t \\ ⩽ {lim inf}_{ε^{'} \to 0} \int_{0}^{T} [| u^{ε^{'}} |^{2} (t) + (C . (X^{ε^{'}} - X_{d}), (X^{ε^{'}} - X_{d})) (t)] d t . \end{matrix} \end{matrix}$ But for any exact control – say $u^{e}$ and $X^{e}$ for the corresponding state function – (i.e. $X^{e} (T) = X_{d} (T)$ ) one has also from the definition of $u^{ε^{'}}$ : $\begin{matrix} (58) & \{\begin{matrix} \int_{0}^{T} [| u^{ε^{'}} |^{2} (t) + (C . (X^{ε^{'}} - X_{d}), (X^{ε^{'}} - X_{d})) (t)] d t \\ ⩽ \int_{0}^{T} [| u^{ε^{'}} |^{2} (t) + (C . (X^{ε^{'}} - X_{d}), (X^{ε^{'}} - X_{d})) (t)] d t \\ + \frac{1}{ε^{'}} (D . (X^{ε^{'}} - X_{d}) (T), (X^{ε^{'}} - X_{d}) (T)) \\ ⩽ \int_{0}^{T} [| u^{e} |^{2} (t) + (C . (X^{e} - X_{d}), (X^{e} - X_{d})) (t)] d t \end{matrix} \end{matrix}$ Therefore, $u^{*}$ is an exact control which minimizes the expression: $\begin{matrix} (59) & u^{e} being an exact control \to \int_{0}^{T} [{| u^{e} |}^{2} (t) + (C . (X^{e} - X_{d}), (X^{e} - X_{d})) (t)] d t \end{matrix}$ where $X^{e}$ is the solution of (1) with the control $u^{e}$ . It is unique because of the strict convexity of the functional defined at (59). Hence, from a standard result, all the sequence $u^{ε}$ converges weakly to $u^{*}$ .

Let us now turn to the strong convergence of $u^{ε}$ to $u^{*}$ . Let us recall that $u^{ε}$ minimizes the criterion (2) and because $u^{*}$ is an exact control, one gets: $\begin{matrix} (60) & \{\begin{matrix} \int_{0}^{T} [| u^{ε} - u^{*} |^{2} + (C . (X^{ε} - X^{*}), (X^{ε} - X^{*}))] \\ = \int_{0}^{T} [| u^{ε} |^{2} + (C . (X^{ε} - X_{d}), (X^{ε} - X_{d}))] \\ - 2 \int_{0}^{T} [u^{ε} u^{*} + (C . (X^{ε} - X_{d}), (X^{*} - X_{d}))] \\ + \int_{0}^{T} {[| u^{*} |^{2} + (C . (X^{*} - X_{d}), (X^{*} - X_{d}))]}^{2} \\ ⩽ 2 \int_{0}^{T} [| u^{*} |^{2} + (C . (X^{*} - X_{d}), (X^{*} - X_{d}))] \\ - u^{ε} u^{*} - (C . (X^{ε} - X_{d}), (X^{*} - X_{d}))], \end{matrix} \end{matrix}$ from the weak convergence result, we deduce that: $\begin{matrix} (61) & \{\begin{matrix} {lim}_{ε \to 0} u^{ε} = u^{*} in the space [L^{2} (] 0, T {[)]}^{m} (p is the number of control) \\ {lim}_{ε \to 0} \sqrt{(C . (X^{ε} - X^{*}), (X^{ε} - X^{*}))} = 0 in the space [L^{2} (] 0, T {[)]}^{N} . \\ But also and more precisely, from the equation satisfied by both X^{ε} and X^{*} : \\ {lim}_{ε \to 0} X^{ε} = X^{*} in the space [H^{1} (] 0, T {[)]}^{N} . \end{matrix} \end{matrix}$ Let us finish this subsection by proving that $u^{0}$ is the exact control in the space $[L^{2} (] 0, T {[)]}^{N}$ which minimizes the strictly convex functional (already introduced at (59)) where $X_{v}$ is the state function associated to the control v: $\begin{matrix} (62) & C (v) = \int_{0}^{T} [| v |^{2} (t) + (C . (X_{v} - X_{d}), (X_{v} - X_{d}))] . \end{matrix}$

Theorem 2.4.
Let us assume that the exact controllability is satisfied. The limit $u^{}$ of the optimal control solution* $u^{ε}$ is the one which minimizes among all the exact control, the quantity: $\begin{matrix} (63) & v \to C (v) = \int_{0}^{T} | v |^{2} + (C . (X_{v} - X_{d}), (X_{v} - X_{d})), \end{matrix}$ where $X_{v}$ is the solution of ( 1 ) with the control v. Therefore, because this solution is unique, one has: $\begin{matrix} u^{0} = u^{}, \end{matrix}$ and this gives an algorithm for computing* $u^{}$ and makes sense to the optimal control formulation (i.e. it is justified if and only if the exact controllability is satisfied otherwise there is no stability of the optimal control* $u^{ε}$ when $ε \to 0$ ).
Proof.
The exact control $u^{0}$ introduced by the asymptotic expansion and characterized previously is such that: $\begin{matrix} (64) & \frac{d X^{0}}{d t} = A X^{0} + F + B u^{0}, X^{0} (0) = X_{0}, X^{0} (T) = X_{d} (T) . \end{matrix}$ Let $U^{e x}$ the non-empty closed convex set of exact controls in $[L^{2} (] 0, T {[)]}^{N}$ . Let us denote by $X_{v}$ the solution of (1) associated to the exact control v. One has: $\begin{matrix} (65) & \{\begin{matrix} \forall v \in U^{e x}, \int_{0}^{T} (u^{0}, v - u^{0}) + (C . (X^{0} - X_{d}), (X_{v} - X^{0})) \\ = - \int_{0}^{T} (^{tr} B . P^{1}, v - u^{0}) + \int_{0}^{T} (C . (X^{0} - X_{d}), (X_{v} - X_{d})) \\ = - \int_{0}^{T} (P^{1}, B (v - u^{0})) + \int_{0}^{T} (C . (X^{0} - X_{d}), (X_{v} - X^{0})) . \end{matrix} \end{matrix}$ We obtain after an integration by parts and because of the equation satisfied by $P^{1}$ and the initial and final conditions satisfied by both $X^{0}$ and $X_{v}$ : $\begin{matrix} (66) & \forall v \in U^{e x}, \int_{0}^{T} (u^{0}, v - u^{0}) + (C . (X^{0} - X_{d}), (X_{v} - X^{0})) = 0 . \end{matrix}$ This is the optimality relation in an affine convex set (as $U^{e x}$ ) which proves that $u^{0}$ is the unique minimizer of (63) (strict convexity of the function: see J. Cea [4]). Finally: $\begin{matrix} u^{0} = u^{*} . \end{matrix}$ □

2.4. A basic solution method for computing $u^{0}$

There are several possibilities in order to compute $u^{0}$ . First of all, one can use a gradient algorithm with optimal step search (see J. Cea [4]), for solving the optimal control problem (i.e. with $ε > 0$ very small but not zero), if we care about minimizing the term in ε in the criterion (2).

If the controllability Assumption 1 is satisfied, then the computational results are correct and the algorithm is not very much parameters sensitive. This is a basic advantage of the exact controllability assumption. Let us also point out that optimal control without exact controllability assumption doesn’t make sense. This is the reason why mathematicians often don’t like the optimal control strategy and they are clearly right. But if the controllability assumption is satisfied, then the algorithms based on the optimal control are fully operational. It remains to the engineers to choose correctly the small parameter ε which can be a tough problem. But there are several possibilities to overcome this difficulty. In the example treated in the next subsection, which is a heating problem of a flat with two rooms but only one radiator, the choice of ε can be performed by comparing the energies spent by the control with the energy required for heating the room.

2.5. A first numerical example with a first order EDO

Let us give a first example using a simple heat equation with $N = 2$ and $m = 1$ (the control is only on the first component of equation (1); i.e. $^{tr} B = [1, 0]$ ). It corresponds to a flat with two rooms exchanging heat but with only one radiator (localized in room 1). The function $X_{d}$ , chosen in this example, aims at trying to give a better heating in one room, for instance the second one, where there is no radiator. and with minimizing the global heating expense (due to the only radiator which is placed in the first room). The temperatures in each room are denoted by $x_{1}$ and $x_{2}$ .

The exact controllability assumption is always satisfied in the following tests. The function $X_{d}$ is represented on the figures by a straight line drawn with dash lines between the two points $x_{1} = x_{2} = 0$ (the initial condition) and $x_{1} = x_{2} = 20$ (the desired target). For each computational test, the walls of the rooms are isotherm (temperature fixed which corresponds – from a mathematical point of view – to homogeneous Dirichlet boundary conditions). But any other conditions could be used. Our goal is just to give an idea of the behavior of the strategy studied in this paper.

The two components of X are plotted on the four sub-figures of each Figs 1–3–5–7 and parametrized by the time. The corresponding limit controls (i.e. for $ε = 0$ ) with respect to the time are plotted for each case on Figs 2–4–6–8. Furthermore, four values of ε have been tested for each case in order to check the numerical convergence. The time T is the same for the three tests and the convergence is accurately obtained for all the tests. The computation of the solutions of the differential equations for X and P has been done using a classical Euler explicit scheme and the stability condition is well satisfied. A θ-Wilson scheme has also been tested but doesn’t lead to improve the results neither concerning the accuracy nor the computational time with larger time steps. Let us turn now to the discussion of the numerical tests.

2.5.1. Test 1: No $X - X_{d}$ term

On the two first Figs 1 and 2 no transient term involving $X - X_{d}$ is taken into account. One can see that the optimal control implies, in this example, an overheating of the first room in order to transfer the energy to the second room.

2.5.2. Test 2: $X - X_{d}$ is applied on $] T / 2, T [$

In the second case, the term $X_{d}$ is used when $t \in] T / 2, T [$ . The results are plotted on Fig. 3 for the trajectory and on Fig. 4 for the control. One can see that the influence of $X_{d}$ is particularly meaningful precisely as expected for $t \in] T / 2, T [$ . But this leads still to a slight over boosting of the temperature in the first room compared to the second one. Hence one can suggest two possibilities: i) increasing the matrix C in order to be closer to the first bisector; ii) restricting the use of $X_{d}$ to a smaller portion of the interval $] 0, T [$ . This is done in the third test because it seems more interesting to reduce the laps of time where $X_{d}$ is influencing the solution. But it would be a special study to determine what is the best choice for both the best time interval.

2.5.3. Test 3: $X - X_{d}$ is applied on $] .8 T, T [$

The results are plotted on Fig. 5 for the temperatures and on Fig. 6 for the exact control. It appears in this case that the energy spent by the control is smaller than in the second case. But the overheating of the first room is more important tan in the previous case. In order to avoid energy spoiling it seems better to start the introduction of the transient term earlier.

2.5.4. Test 4: $X - X_{d}$ is applied on $] 0, T [$

In this last test, the term $X - X_{d}$ is applied on all the interval $] 0, T [$ . Obviously the results are the best. Further more it is worth to notice that the energy spent by the control is not much larger than in the previous cases. But it is larger.

Fig. 1.

No term $X - X_{d}$ . The heating is only applied in the first room. One can see that the control leads to an overheating of the first room which transfer part of its heat to the second room.

Fig. 2.

No $X - X_{d}$ term: the control rise up quickly at the end of the control loop because it is necessary to keep the warming of room 1 and even to over heat it in order to transfer warmth to the second room.

Fig. 3.

A term $X - X_{d}$ is applied between $T / 2$ and T. $X d$ is chosen equal to the straight line between the actual position and the target. The solution is more satisfying because the maximal temperature desired in room 1 is not overtaken. Furthermore the energy expense is not so much increased as one can see by comparing the controls on Figs 2 and 4.

Fig. 4.

The control variable for the simulation case treated on Fig. 3.

Fig. 5.

A term $X - X_{d}$ is applied for $t \in] 0.8 T, T [$ . The heating of the room 2 implies that the transfer of heat is forced from room 1 to room 2 and this implies a decrease of the temperature in room 1. This was not the case if the transient term is applied for $t \in] T / 2, T [$ .

Fig. 6.

The control variable for the simulation case treated on Fig. 5. One can see that it is a little bit too late and the expense of energy is more important than the previous case where the transient term is applied for $t \in] T / 2, T [$ .

Fig. 7.

The term $X - X_{d}$ is applied on all the interval $] 0, T [$ . The heating is only applied in the first room. The results are optimal and the energy spent by the control is similar to the one involved in the previous cases. The temperature increase similarly in the two rooms and the regulation law is simple but not fully obvious.

Fig. 8.

The term $X - X_{d}$ is applied on all the interval $] 0, T [$ . The control rise up regularly up to the end.

2.5.5. Direct computation of the limit of the optimal control for

ε \to 0

Now we discuss hereafter a direct method which is based on the computation of the solution of the equation (30). The idea is to suggest a method which avoids any adjustment on the small parameter ε but also adjustment of the gradient algorithm (choice of the gradient step, initial guess and convergence tests).

Let us first give a matrix representation of the equation (30) using a time step discretization of the equation (1). First of all we consider $N^{t}$ time steps and $N^{t} - 1$ intervals between 0 and T. We set $Δ t = T / (N^{t} - 1)$ . The vectors $X^{0}$ and $δ X^{0}$ are estimated by arrays of dimensions $N \times N^{t}$ (N is the number of components of X). We define the canonical basis of $R^{N}$ by $e_{i}, i = 1, N$ . The component of a vector of $R^{N}$ are denoted by $α^{i}$ . One has for instance: $\begin{matrix} X^{0} = \sum_{i = 1, N} α^{i} e_{i} and δ X^{0} = \sum_{i = 1, N} δ α^{i} e_{i} . \end{matrix}$ But the coefficients $α^{i}$ and $δ α^{i}$ are fonctions of time. At time $j Δ t$ for $j = 0, N^{t} - 1$ we set: $\begin{matrix} α^{i} (j Δ t) = α_{j}^{i}, \end{matrix}$ and the array $α_{j}^{i} \in R^{N} \times R^{N^{t}}$ is denoted by α. The same is done for $δ α$ leading to the array $δ α$ . The vector $Φ \in R^{N}$ represents the initial conditions of $P^{1}$ , and one can write: $\begin{matrix} Φ = \sum_{i = 1, N} Φ^{i} e_{i}, Φ^{i} \in R . \end{matrix}$ Let us denote by $f_{i}$ the functions of time piecewise linear on each segment $[j Δ t, (j + 1) Δ t [$ and such that $f_{i} (j Δ t) = δ_{i j}$ (Kronecker symbol). We also define the basis functions: $\begin{matrix} (67) & F_{i k} (t) = f_{k} (t) e_{i} . \end{matrix}$ Several other definitions are also used. We define by $P (F_{i k})$ or $Q (F_{j k})$ the solutions of: $\begin{matrix} (68) & \{\begin{matrix} \frac{d P}{d t} (F_{i k}) = -^{tr} A . P (F_{i k}) - C . F_{i k}, P (F_{i k}) (0) = 0, \\ \frac{d Q}{d t} (F_{j k}) = -^{tr} A . Q (F_{j k}) - C . F_{j k}, Q (F_{j k}) (0) = 0 . \end{matrix} \end{matrix}$ Due to the linearity of the models involved, one has: $\begin{matrix} (69) & P_{x^{0}}^{1} = \sum_{k = 1, N^{t}} \sum_{i = 1, N} α^{i} (k Δ t) P (F_{i k}) \end{matrix}$ With another respect, if $P (Φ^{i})$ is the solution of: $\begin{matrix} (70) & \frac{d P}{d t} (Φ^{i}) = -^{tr} A . P (Φ^{i}), P (Φ^{i}) = e_{i}, \end{matrix}$ one has: $\begin{matrix} (71) & {\tilde{P}}^{1} = \sum_{i = 1, N} Φ^{i} P (Φ^{i}) . \end{matrix}$

Let us now introduce several matrices in order to define an approximation of problem (30) (we only give the matrices for $N = 2$ and $N_{t} = 4$ for the sake of brevity in the notations). We make use of a reduce integration on $] 0, T [$ for the computation of integrals as follows: $\begin{matrix} (72) & \int_{0}^{T} z (t) d t ≃ \frac{Δ t}{2} [z (0) + z (T) + 2 \sum_{k = 2, N_{t} - 1} z (k Δ t)] . \end{matrix}$ This enables to lump several matrices. Let us set: $\begin{matrix} (73) & G = (\begin{matrix} g & 0 \\ 0 & g \end{matrix}), H = (\begin{matrix} h_{11} & ^{tr} h_{21} \\ h_{21} & h_{22} \end{matrix}) \end{matrix}$ where g is a diagonal matrix with maximum dimensions $N_{t} \times N_{t}$ (hence $4 \times 4$ in our example below) given by (the integral terms involving $X - X_{d}$ are considered for the whole open set $] 0, T [$ but should be limited to the points taken into account in the applications; this restricts significantly the size of the matrix g): $\begin{matrix} (74) & g = Δ t (\begin{matrix} 0.5 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0.5 \end{matrix}) \end{matrix}$ and the sub-matrices $h_{i j}, i, j \in {1, N}$ which components are denoted by $h_{i k}^{k l}$ , are computed by: $\begin{matrix} (75) & h_{i j}^{k l} = \int_{0}^{T} (^{tr} B . P (F_{i k}),^{tr} B . Q (F_{j l})) . \end{matrix}$

Then we also introduce a rectangular matrix – say $R \in R^{N \times N_{t}} \times R^{N}$ – which represents the coupling between Φ and $δ X_{0}$ : $\begin{matrix} (76) & R = {r_{i k}^{j}} = \int_{0}^{T} (^{tr} B . P (F_{i k}),^{tr} B . Q (Φ^{j})), i \in {1, N}, k \in {1, N^{t}}, j \in {1, N} \end{matrix}$ In our example with $N = 2$ and $N^{t} = 4$ the matrix R has the dimensions: 8 lines and 2 columns; ( $(N \times N^{t} = 8) \times (N = 2)$ ). Finally we define the matrix Λ with components $Λ_{k l}, k, l \in {1, N}$ by: $\begin{matrix} (77) & Λ_{k l} = \int_{0}^{T} (^{tr} B . P (Φ^{k}),^{tr} B . Q (Φ^{l})) . \end{matrix}$ The approximation of equation (30) consists in solving: $\begin{matrix} (78) & (\begin{matrix} G + H & R \\ ^{tr} R & Λ \end{matrix}) (\begin{matrix} α \\ Φ \end{matrix}) = S \end{matrix}$ with the last notation for the right hand side $S \in R^{N \times N^{t} + N}$ : $\begin{matrix} (79) & \{\begin{matrix} S = (\begin{matrix} S_{x^{0}} \\ S_{Φ} \end{matrix}), \\ S_{x^{0}}^{i j} = - \int_{0}^{T} (^{tr} B . P_{d},^{tr} B . Q (F_{i j})) - (C . X_{d}, Q (F_{i j})), \\ S_{Φ}^{j} = (X_{0}, e_{j}) - \int_{0}^{T} (^{tr} B . P_{d},^{tr} B . Q (Φ^{j})) . \end{matrix} \end{matrix}$

Remark 4.
In order to decrease the size of the linear system it is worth to restrict the number of point where the integral term involving $X - X_{d}$ is really useful (snapshot points). Once $X_{C}^{0}$ and Φ are computed, one has to compute $P^{1}$ leading to $u^{0}$ and then to $X^{0}$ on all the segment $[0, T]$ .
Remark 5.
In a closed loop setting and for $X_{d}$ equal to zero excepted in the vicinity of $t ≃ T$ , one can use a Schur complement method [10] for solving (78) by eliminating $X_{C}^{0}$ from the first lines of the equation (78). The size of the system to be solved is therefore drastically reduced.
Remark 6.
The numerical results are exactly those obtained with the gradient algorithm for ε very small. Nevertheless, the advantage is that we do not have to adjust the various parameters which appear in the gradient method. The computational time is smaller if $X_{d}$ is non zero on a reduce time interval. But the gradient method is faster if we consider all the interval $] 0, T [$ (which increases the size of the Gammian matrix (78) [17]. Nevertheless. an optimization of the numerical software depending on the computer used would certainly change the conclusion in favor of the direct method.

3. Case of a second order equation

As an example we consider a two dimensional model corresponding to an airfoil which can move in heaving and pitching displacements. In order to fix the ideas, the system is represented on Fig. 9.

3.1. The dynamical model

Fig. 9.

A two degrees of freedom airfoil.

The pitching angle is denoted by α and the heaving by z. A simple linearized model representing the movement of the airfoil can be stated as follows with the notations:

m is the mass of the airfoil;

$J_{0}$ its inertia around point O which is the center of rotation and can only moves in the direction z, $J_{G}$ is the inertia around point G, the center of mass and which is different from O; one has $J_{0} = J_{G} + a^{2} m$ ;

V is the flow velocity far away from the structure;

a is the algebraic distance between O and G;

$c_{z}$ is the lift coefficient given in the Eiffel axis (those of the wind). The corresponding aerodynamical force is $\frac{ϱ S V^{2}}{2} c_{z} e_{z}$ , S being a cross section used as a reference surface;

$c_{m}$ is the pitching coefficient at point O and the pitching moment is $\frac{ϱ S L V^{2}}{2} c_{m}$ , L being a characteristic length.

The equations of this simple model after a linearization around

α = 0

, are:

\begin{matrix} (80) & \{\begin{matrix} m \ddot{z} + a m \ddot{α} + k z = \frac{ϱ S V^{2}}{2} [c_{z} (0) + \frac{\partial c_{z}}{\partial α} (0) α] - m g, \\ J_{0} \ddot{α} + a m \ddot{z} + c α = \frac{ϱ S L V^{2}}{2} [c_{m} (0) + \frac{\partial c_{m}}{\partial α} (0) α] - a m g . \end{matrix} \end{matrix}

In a matrix form one obtains (with $R = \frac{ϱ S V^{2}}{2}$ and $Q = \frac{ϱ S L V^{2}}{2} = R L$ ): $\begin{matrix} \overset{M inertia matrix}{\overset{︷}{(\begin{matrix} m & a m \\ a m & J_{0} \end{matrix})}} (\begin{matrix} \ddot{z} \\ \ddot{α} \end{matrix}) + \overset{A stiffness matrix}{\overset{︷}{(\begin{matrix} k & - R \frac{\partial c_{z}}{\partial α} (0) \\ 0 & c - Q \frac{\partial c_{m}}{\partial α} (0) \end{matrix})}} (\begin{matrix} z \\ α \end{matrix}) = \overset{F forces applied at α = 0}{\overset{︷}{(\begin{matrix} R c_{z} (0) - M g \\ Q c_{m} (0) - M g a \end{matrix})}} \end{matrix}$ The angle $α = 0$ is chosen such that the resulting external forces are equilibrated. Now we add the effect of a flap which is the control (say u) and which acts as both a vertical force and a pitching moment. In fact u is the inclination of the flap and the aerodynamic force induced by this control system is denoted by $B u$ . Therefore it is a one dimensional control function of time even if B is a vector of $R^{2}$ . The dynamical model becomes, with self explanatory notations: $\begin{matrix} (81) & M \ddot{X} + A X = B u, X = (\begin{matrix} z \\ α \end{matrix}), \end{matrix}$ the initial conditions correspond to a perturbation and are written: $\begin{matrix} (82) & X (0) = X_{0}, \dot{X} (0) = X_{1} . \end{matrix}$ It is worth noting that A is not a symmetrical matrix and this is precisely at the origin of a flutter phenomenon (see [7–9]).

3.2. The control problem

Let us first introduce the optimal control criterion function of a control $v \in L^{2} (] 0, T [)$ where X is solution of (81) with the control v and for any $ε > 0$ , by (D is a positive and symmetrical $2 \times 2$ matrix but not necessarily definite): $\begin{matrix} (83) & J^{ε} (v) = \frac{1}{2} [{‖ X (T) ‖}^{2} + {‖ \dot{X} (T) ‖}^{2} + ε \int_{0}^{T} [(D X (s), X (s)) + {| v (s) |}^{2}] d s] . \end{matrix}$ The optimal control model (for $ε > 0$ ) consists in minimizing $J^{ε} (v)$ versus $v \in L^{2} (] 0, T [)$ . Let us point out that in this case the first guess $X_{d}$ is zero but this is just an example and the method could also be applied with a non vanishing first guess as for the first example.

The mathematical analysis is identical to the one given in the previous sections. The main difference is in the solution method which is much easier in this case (because the adjoint state is solution of a stable equation) by solving the optimal control problem with a very small value of the parameter ε. We used an optimal step gradient algorithm. The differential equations are solved again using the central difference scheme and in this case C is arbitrarily chosen equal to the identity (along all the segment $[0, T]$ ). The velocity V of the wind is smaller than the critical ones (there are two). In fact, there are two instabilities: one corresponding to the flutter phenomenon [7–9] and the other one to a negative stiffness of the pitching due to the term: $\begin{matrix} c - Q \frac{\partial c_{m}}{\partial α}, \end{matrix}$ which can be negative if (let us recall that c is the stiffness of the torsion spring): $\begin{matrix} V ⩾ \sqrt{\frac{2 c}{ϱ S L \frac{\partial c_{m}}{\partial α}}} . \end{matrix}$ On Fig. 10 we have plotted the solution of the exact control with $C = 0$ (we used a direct solver as explained in previous sections). The control and its derivative versus the time have also been plotted. On Fig. 11 is represented the exact control (computed with ε very small and a gradient algorithm) with $C = I_{d}$ . One can see that the control is still exact as scheduled by the theory, but the magnitude of the variables z and α are smaller than in the previous case. Obviously the magnitude of the control is larger but this is not a difficulty as far as it only concerns the amplitude of the pitching movement of the rear flap. In both cases the scales used on the pictures are arbitrary and the aerodynamic coefficients used are those of a classical airfoil obtained from our wind tunnel similar to those of a NACA-0012 airfoil [1].

Fig. 10.

Top left and right are respectively the heaving and the angle of attack of the airfoil. At the bottom left is plotted the control (oscillation of the rear flap) and its derivative versus time on the right (T is in 0.001 s, α in $0.1 rad$ and z in $0.01 Meter$ . The wind velocity is $20 m/s$ ).

Fig. 11.

The solution of the limit optimal control model ( $ε = 0$ ) is plotted on these four pictures. One can see that the magnitude of the heaving (top left) and of the angle of attack are smaller. Nevertheless the control is still exact as announced in the theoretical analysis. The new control is plotted on the bottom left picture and it is larger than in the previous case ( $C = 0$ ). Because we used the optimal control method the decay of the control criterion dring the gradient algorithm has been plotted on the bottom right picture.

4. Conclusion

In this paper we have discussed the asymptotic behavior of an optimal control with a quadratic criterion involving both the final values and a transient term of a state function, when the small parameter ε (representing the marginal cost of the control) tends to zero. The new point is that it is possible to introduce in the cost of the control any quadratic and positive terms which can take into account the values of a state function during the control process. For instance, one can introduce intermediate targets for the state function. When the marginal cost parameter ε tends to zero and if the exact controllability of the system is satisfied, the optimal control converges to the cheapest exact control in term of the transient term in the cost function which has been defined. Furthermore, there is a way to compute directly this exact control avoiding any choice of ε. But in some case, it is more convenient to use the optimal control algorithm (i.e. $ε > 0$ ). Two examples have been discussed one for a first order ordinary differential equation (control of the heating of a flat with two rooms and only one radiator) and the second one is a stabilization of an airfoil which has two degrees of freedom. In both cases the influence of the additional transient term is clearly visible. In fact it is worth to underline that this strategy is an extension to ordinary differential equations of well known model used by automation engineers and named phase control [18]. This last one is a finite dimensional formulation of the so-called HUM algorithm developed by J.L. Lions [16] for partial differential equations.

References

http://airfoiltools.com/airfoil/details?airfoil=n0012-il.

Bellman, Dynamic Programming, 1956, reedited by Dover Books on Computer Science, edited in 2003.

Bellman, Introduction to the Mathematical Theory of Control Process, Vol. 1, Academic Press, Cambridge MA, USA, 1958.

Cea, Optimisation: Théorie et Algorithmes, Dunod, Paris, 1968.

Descombes, Integration, Hermann, Paris, 1973.

Destuynder, Analyse et Contrôle des Équations Différentielles, Hermès-Lavoisier, Paris-Londres, 2010.

Destuynder and

Fabre, An introduction to quasi-static aeroelasticity, in: Numerical Simulation and Engineering Trends and Applications, SEMA-SMAI Lecture Notes of the XVIII Jacques-Louis Lions Spanish-French School, Springer Series, Vol. 24, Berlin, 2021.

Faurre and

Robin, Eléments d’Automatique, Dunod, Paris, 1984.

Y.C.

Fung, Introduction to the theory of aeroelasticity, in: Dover Books on Aeronautical Engineering, Dover Editions, New York, 1969.

10.

R.A.

Horn and

C.R.

Johnson, Matrix Analysis, Cambridge University Press, New York, 2012.

11.

R.E.

Kalman, Contributions to the theory of optimal control, Bol. Soc. Mat. Mexicana 5 (1960), 102–119.

12.

R.E.

Kalman, On the general theory of control, in: Proc. of the First Int. Cong. on Auto. Control, Vol. 1, 2nd edn, Butterworth Scientific Inc., London, 1964.

13.

R.E.

Kalman, When is a linear control system optimal?, J. Basic Eng. 86(1) (1964), 51–60. doi:10.1115/1.3653115.

14.

J.L.

Lions, Contrôle Optimal de Systèmes Gouvernés Par des Équations aux Dérivées Partielles, Dunod, Paris, 1969.

15.

J.L.

Lions, Perturbations Singulières dans les Problèmes aux Limites et en Contrôle Optimal, Lectures Notes in Mathematics, Vol. 323, Springer, Berlin, Heidelberg, 1973.

16.

J.L.

Lions, Contrôlabilité Exacte, Perturbation et Stabilisation de Systèmes Distribués, RMA, Masson, Paris, 1988.

17.

Ogata, Modern Control Engineering, 4th edn, Prentice Hall International Editions, New York, 2001.

18.

Pallu de la Barrière, Cours d’Automatique Théorique, Dunod, Paris, 1966.

19.

L.S.

Pontryagin,

V.G.

Boltyanskii,

R.V.

Gamkrelidze and

E.F.

Mishchenko, The Mathematical Theory of Optimal Processes, Interscience, 1962.

20.

A.N.

Tikhonov, On the stability of inverse problems, Dokl. Akad. Nauk SSSR 39(5) (1943), 195–198.

21.

A.N.

Tikhonov and

V.A.

Arsenin, Solution of Ill-Posed Problems, Winston and Sons, Washington, 1977 (ISBN 0-470-99124-0).

Asymptotic method and transient terms in exact controls

Abstract

Keywords

1. Introduction

2. Characterization of u ε

Assumption 1 (Controllability (see R. Bellman [2] or I. Pontryagin [19])).

1 ∀ a , b ∈ R , α ∈ R + : 2 a b ⩽ α a 2 + 1 α b 2 .

2.5. A first numerical example with a first order EDO

2.5.1. Test 1: No X − X d term

2.5.2. Test 2: X − X d is applied on ] T / 2 , T [

2.5.3. Test 3: X − X d is applied on ] .8 T , T [

2.5.4. Test 4: X − X d is applied on ] 0 , T [

3.1. The dynamical model

References

2. Characterization of $u^{ε}$

¹
$\forall a, b \in R$ , $α \in R^{+} : 2 a b ⩽ α a^{2} + \frac{1}{α} b^{2}$ .

2.5.1. Test 1: No $X - X_{d}$ term

2.5.2. Test 2: $X - X_{d}$ is applied on $] T / 2, T [$

2.5.3. Test 3: $X - X_{d}$ is applied on $] .8 T, T [$

2.5.4. Test 4: $X - X_{d}$ is applied on $] 0, T [$