<figure><img src="image_1.png" /></figure>

HAL
open science

# Indirect adaptive control of unknown diffusion equation

Tudor-Bogdan Airimitoaie, Christophe Collewet

► To cite this version:

Tudor-Bogdan Airimitoaie, Christophe Collewet. Indirect adaptive control of unknown diffusion equation. 7th AIAA Flow Control Conference, Jun 2014, Atlanta, Antarctica. pp.9. hal-01101083

HAL Id: hal-01101083

https://inria.hal.science/hal-01101083

Submitted on 30 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Indirect adaptive control of unknown diffusion equation

Tudor-Bogdan Airimițoaie and Christophe Collewet

IRSTEA / INRIA Rennes - Bretagne Atlantique, Fluminance team, Rennes, France.

This paper concerns the control of a diffusion equation with unknown or uncertain diffusivity. The purpose is to demonstrate the capabilities of adaptive control based on the identification of the unknown diffusion coefficient. This problem is similar to the control of flow with unknown Reynolds number which will be investigated in a future publication. The present version gives only the main ideas which will be fully detailed in the final version of the paper.

## Notation

We use the subscript notation for partial derivative, for example $\xi_x(x, y)$ represents the derivative of the dependent variable $\xi(x, y)$ with respect to the independent variable $x$, $\partial\xi/\partial x$, and similarly $\xi_{xy}(x, y) = \partial^2\xi/\partial x\partial y$. $\xi_{i,j}$ denotes the value of a dependent variable in a point $(i\Delta x, j\Delta y)$ on a grid of discretized positions, where $\Delta x$ and $\Delta y$ represent the discretization steps in each of the two independent variables.

## I. Introduction

In the control of systems from all areas of engineering, model based designs are often found. These are based on computing a control law from knowledge upon the system given by some kind of modeling of the input-output behavior. However, in practical control applications, there often exist unknown or uncertain parameters. This is the case also in the problem of flow control where various parameters are unknown, depending on the characteristics of the system, and can change their values over time. Therefore, robust or adaptive control techniques have been used to improve performances with respect to uncertainties in the model of the system.

In model based control (MBC), uncertainties in the modeling of the system play an important role on performances. They can influence the stability of the closed loop system, the reference tracking capabilities or disturbance rejection efficiency. Two main approaches have been developed for dealing with uncertainties. One of them is robust control and it is based on designing a linear controller capable of maintaining good performances under reasonable parametric uncertainties. The second is adaptive control,¹ which tries to estimate the uncertainty and update the controller in order to adapt to the changing environment. The first approach is popular for its simplicity, but the second is more efficient when the level of uncertainty is important, as is the case in this article. In fact, we study the reference tracking capability under an uncertain parameter (the unknown diffusivity) and, as it will be shown in the simulation section, the static gain matrix plays an important role here. However, when the diffusivity parameter changes its value, the gain matrix needs to be updated, otherwise very large differences from the desired values are obtained. For this purpose, it is necessary to estimate the true value of the diffusivity and update the controller.

The presented method is a step towards the adaptive control of flow with unknown or time-varying Reynolds number. But different from the usual objective of stabilization around an equilibrium position, as is the case in various flow control papers<sup>2,3,4,5,6</sup>, we are interested here in tracking a certain reference signal<sup>7</sup>. It is well known that the robustness of flow control algorithms is very dependent on this parameter and instability can occur if the linearized system is "too far away", from a robust stability point of view, from the real system<sup>8</sup>. A solution can be the design of adaptive or gain scheduled controllers but some of the parameters of the true system have to be estimated. The solution proposed in this paper for the estimation of a diffusion equation provides a promising solution for this problem.
In the next sections of this article, we describe the classical approach to modeling infinite dimensional systems, as it is the case in flow control, which is spatial discretization using a carefully chosen set of grid points. A basic method of approximation is presented in section II and applied in section III to the discretization of the diffusion equation with time-varying diffusivity. Then, in section IV, we describe the indirect adaptive control approach that is the main idea presented in this article. Some simulation results are shown in section V. Finally, section VI concludes the paper.

## II. Finite-difference approximation method

Most of the systems encountered in nature can be described by partial differential equations (PDEs). Therefore, in an attempt to control such systems, the need to better understand them arises. As such, numerical analysis has seen a lot of development in the last decades as the computational power of personal computers increased. Numerical analysis is based upon discretisation of PDEs and various techniques have been described.

In this section, a brief review of the finite-difference method (FDM) is presented (for a more in depth description see⁹). This method is based on approximating derivatives by differences, resulting in an algebraic representation of the PDEs. The first step is to replace the continuous domain by a finite-difference mesh. For a two-dimensional (2-D) problem, $\xi(x, y)$ is replaced by $\xi(i\Delta x, j\Delta y)$ (written also $\xi_{i,j}$) for given distance increments $\Delta x$ and $\Delta y$.

Recall next the definition of the derivative of a function $\xi(x, y)$ at $x = x_0, y = y_0$:

$$ \frac{\partial \xi}{\partial x} = \lim_{\Delta x \to 0} \frac{\xi(x_0 + \Delta x, y_0) - \xi(x_0, y_0)}{\Delta x} \quad (1) $$

Another way to find this equation is by considering first a Taylor-series expansion for $\xi(x_0 + \Delta x, y_0)$ about $\xi(x_0, y_0)$ which is given by

$$ \xi(x_0 + \Delta x, y_0) = \xi(x_0, y_0) + \left. \frac{\partial \xi}{\partial x} \right|_0 \Delta x + \left. \frac{\partial^2 \xi}{\partial x^2} \right|_0 \frac{(\Delta x)^2}{2!} + \dots \quad (2) $$

or in the $(i, j)$ notation

$$ \xi_{i+1,j} = \xi_{i,j} + \left. \frac{\partial \xi}{\partial x} \right|_{i,j} \Delta x + \left. \frac{\partial^2 \xi}{\partial x^2} \right|_{i,j} \frac{(\Delta x)^2}{2!} + \dots \quad (3) $$

which gives

$$ \left. \frac{\partial \xi}{\partial x} \right|_{i,j} = \frac{\xi_{i+1,j} - \xi_{i,j}}{\Delta x} + O(\Delta x), \quad (4) $$

where $(\xi_{i+1,j} - \xi_{i,j})/\Delta x$ represents the "forward" finite-difference approximation of the derivative and $O(\Delta x)$ is the truncation error between the derivative and the proposed representation.

In a similar manner, it is also possible to find a "backward" formula

$$ \left. \frac{\partial \xi}{\partial x} \right|_{i,j} = \frac{\xi_{i,j} - \xi_{i-1,j}}{\Delta x} + O(\Delta x) \quad (5) $$

and a "central" formula

$$ \left. \frac{\partial \xi}{\partial x} \right|_{i,j} = \frac{\xi_{i+1,j} - \xi_{i-1,j}}{2\Delta x} + O(\Delta x)^2. \quad (6) $$

As for the second derivative, this can be evaluated using the central difference formula

$$ \left. \frac{\partial^2 \xi}{\partial x^2} \right|_{i,j} = \frac{\xi_{i+1,j} - 2\xi_{i,j} + \xi_{i-1,j}}{(\Delta x)^2} + O(\Delta x)^2. \quad (7) $$

## III. Diffusion equation

Consider the following diffusion equation, described by the initial boundary value problem (IBVP):

$$ \xi_t(x, t) = \theta^* \xi_{xx}(x, t) + b^T(x)u(t), \quad x \in (0, 1), \ t > 0, \quad (8) $$
with boundary conditions (BCs)

$$ \xi_x(0, t) = \xi_x(1, t) = 0, \qquad (9) $$

and initial condition (IC)

$$ \xi(x, 0) = \xi_0(x) \qquad (10) $$

where $\xi(x, t)$ represents the system state, $u(t)$ is a control signal, and $\theta^*$ is a positive unknown and possibly time-varying parameter that represents the diffusivity. This is also known as the heat equation with time-varying thermal conductivity.

The input influences the PDE (8) only at certain locations specified by the input shaping function

$$ b^T(x) = [b_1(x) \ b_2(x) \ \dots \ b_{10}(x)] \qquad (11) $$

with

$$ b_{i+1}(x) = \begin{cases} 1, & 0.1i < x < 0.1 + 0.1i \\ 0, & \text{otherwise,} \end{cases} \qquad (12) $$

$$ \forall i \in \{0, 1, 2, 3\}. $$

The output is described by the equation

$$ y(t) = \int_{0}^{1} \xi(x, t)c(x)dx, \quad t > 0. \qquad (13) $$

where $c(x)$ is given by

$$ c^T(x) = [c_1(x) \ c_2(x) \ \dots \ c_{80}(x)] \qquad (14) $$

with

$$ c_{i+1}(x) = \begin{cases} 1, & \frac{i}{80} < x < \frac{i+1}{80} \\ 0, & \text{otherwise} \end{cases} \qquad (15) $$

The eigenvalues of the PDE (8) can be computed as

$$ \lambda_i = -\theta^*(i-1)^2\pi^2, \quad i \ge 1. \qquad (16) $$

Note that the system has one zero eigenvalue and an infinity of negative ones. The numerical simulation of the IBVP can be done by using the method of lines. The semi-discretisation of the PDE is done on the spatial domain by dividing the interval on which $x$ is defined into $N-1$ equal subintervals, obtaining thus $N$ grid points. The method of lines is used in initial boundary value problems and it consists in discretising the system along all minus one of the variables. This results in a system of ordinary differential equations (ODE) which can be solved using various existing methods (see for example the Euler method and the various alternatives of the Runge-Kutta method described in¹⁰) efficiently implemented in simulation software such as MATLAB, Maple, Mathematica, or Scilab.

After discretization and taking into account the boundary conditions, a lumped parameter system described by a state space representation of the form

$$ \dot{\chi}(t) = \theta^* A_L \chi(t) + B_L u(t), \qquad (17a) $$

$$ y_L(t) = C_L \chi(t), \qquad (17b) $$

is obtained, where $\chi \in \mathbb{R}^N$, $A_L \in \mathbb{R}^{N,N}$

$$ A_L = \frac{1}{\Delta x^2} \begin{pmatrix} -2 & 2 & 0 & 0 & 0 & \dots & 0 \\ 1 & -2 & 1 & 0 & 0 & \dots & 0 \\ 0 & 1 & -2 & 1 & 0 & \dots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \dots & 0 & 1 & -2 & 1 & 0 \\ 0 & \dots & 0 & 0 & 1 & -2 & 1 \\ 0 & \dots & 0 & 0 & 0 & 2 & -2 \end{pmatrix}, \qquad (18) $$
$B_L \in \mathbb{R}^{N,10}$, and $C_L \in \mathbb{R}^{N,N}$ have zero values everywhere exception of a few elements of value 1 at grid positions corresponding to eqs. (12) and (15). Note that with the given output shaping function, $C_L = I_N$ and therefore $y_L(t) = \chi(t)$. The diffusivity is considered as time-varying and the spatial discretization step is $\Delta x = \frac{1}{N-1}$. Note that, even though $y(t)$ and $y_L(t)$ have the same dimensions, they differ as $y(t)$ represents the true output from the distributed parameters systems while $y_L(t)$ is the output from the space discretized system and represents an approximation of $y(t)$.

## IV. Estimation and adaptive control

We are interested here in the adaptive control of the unknown or uncertain diffusion equation in order to follow a prespecified profile. The following control is considered

$$u(t) = u_r(t) + u_{\chi}(t) = G \cdot r(t) - K \cdot y(t), \quad (19)$$

where $G$ is designed to assure unit static gain transfer from input $r(t)$ to output $y(t)$, and $K$ is a state feedback gain, as shown also in figure 1. In the next subsections we further detail the controller design and identification methods.

### A. Linear controller design for reference tracking

From the state space representation it is straightforward to obtain the input-output transfer matrix of the closed loop system with feedback gain $K$ and feedforward gain $G^{11}$

$$H_{CL}(s) = G \cdot C(s \cdot I_N - \theta^* A_L + B_L \cdot K)^{-1} B \quad (20)$$

and using the final value theorem, it is easy to compute $G$ for unit static gain from input to output as

$$G = -(C(\theta^* A_L - B_L \cdot K)^{-1} B)^{-1}. \quad (21)$$

<figure><img src="image_2.png" /><figcaption>Figura 1. Closed loop state feedback control scheme with unit static gain.</figcaption></figure>

Regarding the state feedback gain matrix $K$, a linear quadratic regulator (LQR) design method$^{12}$ is used in this paper. More precisely, it is an optimal controller design methodology for continuous time state space systems of the form (17), which allows to design a control of the form

$$u_{\chi}(t) = -K\chi(t) \quad (22)$$

that minimizes a quadratic cost function of the form

$$J = \int_{0}^{\infty} (\chi^{T} Q \chi + u^{T} R u) dt. \quad (23)$$

In the last equation, $Q$ and $R$ are state and input weighting matrices chosen in order to satisfy certain design restrictions. The solution of the state feedback gain matrix is given by

$$K = R^{-1} B_L^T P, \quad (24)$$

where $P$ is found by solving a continuous time algebraic Riccati equation

$$A_L^T P + P A_L - P B_L R^{-1} B_L^T P + Q = 0. \quad (25)$$
## B. Identification of the unknown diffusivity

As already mentioned in the introduction, the estimation of the unknown parameter plays an important role and can improve performances by use in an adaptive control scheme. In this subsection, we detail the algorithm used for the estimation of the unknown diffusivity. The algorithm development is done under the assumption of unknown but positive and constant diffusivity $\theta^*$ (as it is also the case for the Reynolds number)

$$ \dot{\chi}(t) = \theta^* A_L \chi(t) + B_L u(t). \qquad (26) $$

The idea is that of considering an estimator of the unknown system (where the unknown is represented by the parameter $\theta^*$)

$$ \dot{\chi}(t) = \hat{\theta}(t) A_L \hat{\chi}(t) + B_L u(t), \qquad (27) $$

$\hat{\theta}(t)$ denoting the estimated value of $\theta^*$ at time $t$. The difference between the true and the estimated parameter defines the estimation error

$$ \tilde{\theta}(t) = \theta^* - \hat{\theta}(t). \qquad (28) $$

The state estimation error is then given by

$$ \epsilon(t) = \chi(t) - \hat{\chi}(t) \qquad (29) $$

and its time derivative can be calculated as

$$ \dot{\epsilon}(t) = \dot{\chi}(t) - \hat{\chi}(t) = \theta^* A_L \chi(t) - \hat{\theta}(t) A_L \hat{\chi}(t). \qquad (30) $$

A Lyapunov design approach will be used next to estimate the unknown parameter. This technique is often used in the analysis of the stability of dynamic systems. The basic idea consists in finding a Lyapunov function which represents the energy of the system. This function has to be positive definite everywhere except at the origin. If furthermore its time derivative is negative definite, the origin of the system is asymptotically stable. In our case, the origin is defined as the unknown parameter $\theta^*$, hence estimation of this parameter can be realized by designing an appropriate Lyapunov function and update law for $\hat{\theta}(t)$. When the time derivative of the Lyapunov function is not negative definite, a persistence of excitation condition¹³ has to be added, which basically means that the input of the system is reach enough in frequency in order to obtain at the output the necessary information needed for estimation.

Consider a Lyapunov candidate function

$$ V = \epsilon(t)^T P \epsilon(t) + \tilde{\theta}^2(t), \qquad (31) $$

where $P$ is chosen such that

$$ A_L^T P + P A_L = -Q \qquad (32) $$

and $Q$ is any $N \times N$ symmetric positive definite matrix.¹ One can see that $V$ is positive semi-definite, because $P$ obtained from (32) is positive definite. A persistent excitation condition has to be considered, as will be detailed next, to ensure that $V$ is equal to zero only at the origin ($\hat{\theta}(t) = \theta^*$).

The time derivative of $V$ is then given by

$$ \dot{V} = \epsilon^T(t) P \dot{\epsilon}(t) + \dot{\epsilon}^T(t) P \epsilon(t) + 2\tilde{\theta}(t) \dot{\tilde{\theta}}(t) \qquad (33) $$

and taking into account that $\dot{\epsilon}(t)$ can be manipulated by adding and subtracting $\theta^* A_L \hat{\chi}(t)$ into

$$ \dot{\epsilon}(t) = \theta^* A_L \chi(t) - \hat{\theta}(t) A_L \hat{\chi}(t) + \theta^* A_L \hat{\chi}(t) - \theta^* A_L \hat{\chi}(t) = \theta^* A_L \epsilon(t) + \tilde{\theta}(t) A_L \hat{\chi}(t). \qquad (34) $$

then $\dot{V}$ becomes

$$ \dot{V} = \theta^* \epsilon^T(t) (A_L^T P + P A_L) \epsilon(t) + 2\tilde{\theta}(t) \epsilon^T(t) P A_L \chi(t) + 2\tilde{\theta}(t) \dot{\tilde{\theta}}(t). \qquad (35) $$

From the last equation, in order to eliminate the last two terms, it is useful to consider the update law for the estimation of $\theta^*$ as

$$ \dot{\tilde{\theta}}(t) = -\epsilon^T(t) P A_L \chi(t), \qquad (36) $$

leading to

$$ \dot{\tilde{\theta}}(t) = \frac{d}{dt} (\theta^* - \hat{\theta}(t)) = -\dot{\hat{\theta}}(t) = \epsilon^T(t) P A_L \chi(t). \qquad (37) $$
Taking into account (35), (32), and the fact that $\theta^*$ is always positive, it remains that $\dot{V} \le 0$. With this result, the Lyapunov stability theorem can be applied taking into account also a persistent excitation (PE) condition to conclude on the convergence of the proposed estimator¹. This is done in order to complement the fact that the time derivative of the Lyapunov function is not negative definite but negative semi-definite. Usually in adaptive control and parametric identification, the PE condition is used in order to ensure that the convergence to zero of the estimation error is not due to a some unfortunate combination of input signals but is due to the convergence of the parameters to the true ones. In our case, for the estimation of a single parameter, the PE condition resumes to the input being different than zero. Therefore, a step signal is sufficient for estimation and this condition is always satisfied.

**Remark 1:** this estimator has been developed under the assumption of known system matrices $A_L$ and $B_L$. Although in reality these are not known exactly, sufficiently good approximations can be obtained which will ensure that $\dot{V} \le 0$.

**Remark 2:** current work concerns the design of a linear parameter varying (LPV) control law with the diffusivity as time-varying parameter. The most important difficulty being that in LPV control measurement instead of estimation is assumed for the scheduling parameter. This would however play nicely in flow control, where, given for example the control of the Poiseuille flow, it would allow design of LPV controllers using as time-varying scheduling parameter the Reynolds number.

### C. Indirect adaptive control

For the purpose of tracking a certain reference signal, $r(t)$, it is important to compute update the matrix $G$ in real time so that unit static gain from input to output is assured, otherwise the reference signal can not be tracked efficiently.

Using the described parametric estimation procedure, we can obtain the value of the diffusivity and update the gain matrix from (21)

$$ \hat{G}(t) = - \left( C \left( \hat{\theta}(t) A_L - B_L \cdot K \right)^{-1} B \right)^{-1}. \quad (38) $$

### V. Simulation results

In this section we consider the control of temperature on a one dimensional heat conductor which is one of the possible systems that can be described through the IBVP (8)-(10). We consider the non-dimensional equation and a spatial discretisation with $N = 120$ (further reduction of this model can be obtained using standard model reduction techniques from the literature). To simulate the true system, we consider a finite difference discretization with $N = 600$ points.

<figure><img src="image_3.png" /><figcaption>Figura 2. Comparison of desired and obtained temperature profiles for the unknown diffusion system (non-adaptive case).</figcaption></figure>

In the following simulation examples we motivate the interest of using the proposed adaptive control technique. First, in figure 2, the obtained temperature profile is shown in the non-adaptive case. The large
difference between the desired and the measured profiles comes from the fact that the gain matrix $G$ is computed for a 0.1 value of the diffusivity, while the true one is 3.

The behavior in a realistic setup, where the real process has been simulated using a finite-difference discretization with $N = 600$ points, while for the controller design the same model as before (with $N = 120$) has been used. In this case, the matrices $A_L, B_L$, and $C_L$ represent approximations of the true system. The diffusivity estimation result is shown in figure 4. A small deviation from the true value is observed which has as a consequence a deviation from the desired reference as can be seen in figure 3.

<figure><img src="image_4.png" /><figcaption>Figura 3. Comparison of desired and obtained temperature profiles for the unknown diffusion system - realistic setup (adaptive case).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figura 4. Estimation of the diffusivity - realistic setup (the true value is 3, the initial value is 0.1).</figcaption></figure>

For the adaptive case and perfect modeling (that is, the model and the true system are represented by the same $A_L, B_L$, and $C_L$), one obtains the profile as shown in figure 5. The reference signal is tracked perfectly at the output of the system. This is because the diffusivity is estimated very well, see figure 6, allowing thus to update the gain matrix $G$ to its correct value.

## VI. Concluding remarks

In this article, we have presented an adaptive control approach for reference tracking on an unknown diffusivity equation based on a first step of estimating the unknown parameter and then updating a static gain matrix. A Lyapunov design approach is used and convergence of the parameter estimator is shown. Simulation results are given to validate the design. It is believed that this method can be further developed for use in control of flow with unknown Reynolds number.
<figure><img src="image_6.png" /><figcaption>Figura 5. Comparison of desired and obtained temperature profiles for the unknown diffusion system (adaptive case).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figura 6. Estimation of the diffusivity (the true value is 3, the initial value is 0.1).</figcaption></figure>
Bibliografie

¹Ioannou, P. A. and Sun, J., *Robust adaptive control*, Dover Publications, 2012.

²Ravindran, S., "Adaptive reduced-order controllers for a thermal flow system using proper orthogonal decomposition," *SIAM Journal on Scientific Computing*, Vol. 23, No. 6, 2002, pp. 1924–1942.

³Vazquez, R. and Krstic, M., "A closed-form feedback controller for stabilization of the linearized 2-D Navier-Stokes Poiseuille system," *Automatic Control, IEEE Transactions on*, Vol. 52, No. 12, 2007, pp. 2298–2312.

⁴Aamo, O. and Fossen, T., "Tutorial on Feedback Control of Flows, Part I: Stabilization of Fluid Flows in Channels and Pipes," *Modeling, Identification and Control*, Vol. 23, No. 3, 2002, pp. 161–226.

⁵Joshi, S., Speyer, J., and Kim, J., "A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow," *Journal of Fluid Mechanics*, Vol. 332, 1997, pp. 157–184.

⁶McKernan, J., Papadakis, G., and Whidborne, J. F., "A Linear State-Space Representation of Plane Poiseuille Flow for Control Design- A tutorial," *International Journal of Modelling, Identification and Control*, Vol. 1, No. 4, 2006, pp. 272–280.

⁷Cochran, J. and Krstic, M., "Motion planning and trajectory tracking for three-dimensional Poiseuille flow," *Journal of Fluid Mechanics*, Vol. 626, 2009, pp. 307–332.

⁸Kim, J. and Bewley, T., "A linear systems approach to flow control," *Annual Review of Fluid Mechanic*, Vol. 39, No. 1, January 2007, pp. 383–417.

⁹Anderson, D. A., Tannehill, J. C., and Pletcher, R. H., *Computational Fluid Mechanics and Heat Transfer*, McGraw-Hill, 1984.

¹⁰Butcher, J. C., *Numerical Methods for Ordinary Differential Equations, Second Edition*, John Wiley & Sons Ltd, 2008.

¹¹Skogestad, S. and Postlethwaite, I., *Multivariable feedback control: analysis and design*, John Wiley & Sons Ltd, 2005.

¹²Zhou, K., Doyle, J. C., Glover, K., et al., *Robust and optimal control*, Vol. 40, Prentice Hall New Jersey, 1996.

¹³Ljung, L., "System identification: theory for the user," *Prentice Hall information and system sciences series*, 1999.