Data-driven reference model selection and application to L-DDC design

TL;DR

This paper introduces a data-driven method for selecting reference models in control design, ensuring they are both achievable by the plant and meet desired performance criteria, specifically for LTI systems.

Contribution

It proposes a novel approach to build reproducible and performance-oriented reference models using plant instability estimation and data-driven stability analysis.

Findings

Improved reference model selection enhances control performance.

The method enables validation via small-gain theorem.

Application to L-DDC demonstrates practical effectiveness.

Abstract

The choice of a reference model in data-driven control techniques is a critical step. Indeed, it should represent the desired closed-loop performances and be achievable by the plant at the same time. In this paper, we propose a method to build such a reference model, both reproducible by the system and having a desired behaviour. It is applicable to Linear Time-Invariant (LTI) monovariable systems and relies on the estimation of the plant's instabilities through a data-driven stability analysis technique. The L-DDC (Loewner Data Driven Control) algorithm is used to illustrate the impact of the choice of the reference model on the control design process. Finally, the proposed choice of specifications allows to use a controller validation technique based on the small-gain theorem.

Tables

Table 1. Table 1: Estimation of the RHP poles of the plant.

RHP poles of $P_{502}$	$3.83\times {10}^{-5}\pm 0.848\times {10}^{-2}ı$
Estimated RHP poles	$1.07\times {10}^{-4}\pm 0.852\times {10}^{-2}ı$
Estimated RHP poles in [13]	$0.99\times {10}^{-4}\pm 0.89\times {10}^{-2}ı$

Equations

$$\left\{\begin{array}[]{l}\forall i=1\dots n_{z},\ M(z_{i})=0\\ \forall j=1\dots n_{p},\ M(p_{j})=1\end{array}\right.,$$

$$B_z(s)=\prod _{i=1}^{n_z}\frac{s-z_i}{s+z_i}$$

$$B_p(s)=\prod _{j=1}^{n_p}\frac{s-p_j}{s+p_j}$$

$$F(s)=\frac{{\sum }_{k=1}^{n_p}{\gamma }_k{l}_k(s)}{{\prod }_{j=1}^{n_p}(s+p_j)},$$

$$P(ı\omega )=P^s(ı\omega )+P^{as}(ı\omega )$$

$$M(s)=\frac{1}{1+\tau s},\tau =1s.$$

$$K_2(s)=\frac{39.084(s^2+0.04163s+0.003132)}{s(s+0.002751)}$$

$$K_{[14]}(s)=\frac{54.47s^2+2.317s+0.02446}{s^2+0.002033s+4.374e-6}.$$

$$K_{M_{14}}(s)=\frac{28.282(s+0.08652)(s+0.00982)}{(s^2+0.001967s+4.192e-06)}.$$

$$\begin{array}[]{ccl}h(x,s,Q_{0})&=&G_{e}(x,s,Q_{0})q_{e}(s)+G_{s}(x,s,Q_{0})q_{s}(s)\\ &=&P(x,s,Q_{0})\left[\begin{array}[]{c}q_{e}\\ q_{s}\end{array}\right]\end{array},$$

$$\begin{array}[]{r@{=}l}G_{e}(x,s,Q_{0})&\frac{\lambda_{1}(s)e^{\lambda_{2}(s)L+\lambda_{1}(s)x}-\lambda_{2}(s)e^{\lambda_{1}(s)L+\lambda_{2}(s)x}}{B_{0}s(e^{\lambda_{1}(s)L}-e^{\lambda_{2}(s)L})}\\ G_{s}(x,s,Q_{0})&\frac{\lambda_{1}(s)e^{\lambda_{1}(s)x}-\lambda_{2}(s)e^{\lambda_{2}(s)x}}{B_{0}s(e^{\lambda_{1}(s)L}-e^{\lambda_{2}(s)L})}\end{array}$$

$$M(s)=\frac{1}{1+\frac{2\xi }{{\omega }_0}s+\frac{s^2}{{\omega }_0^2}},$$

$$\forall i=1\dots N,K(ı{\omega }_i)=\frac{M(ı{\omega }_i)}{{\Phi }_i-M(ı{\omega }_i){\Phi }_i}.$$

$$K(s)=\frac{3.3492e-05(s^2-0.0663s+0.007729)}{(s+2.001e-05)(s+0.001161)}.$$

Full text

Appendix to [1]:

Data-driven reference model selection

and application to L-DDC design

Pauline Kergus1⋆, Martine Olivi2, Charles Poussot-Vassal1 and Fabrice Demourant1 Thanks to the whole FACTAS team at INRIA and particularly to Adam Cooman for his precious advice on data-driven stability analysis.1ONERA, DTIS, Information Processing and Systems, Toulouse, France2FACTAS team, INRIA Sophia Antipolis, France⋆Corresponding author: pauline.kergus@onera.fr

1 Introduction

For many applications, a mathematical description of the system, derived from physical laws, is not available. In this case, the controller has to be designed on the basis of experimental measurements. The first solution consists in identifying a model of the plant and then using any kind of model-based technique to obtain a control law (indirect methods). It is indicated for problems where a reliable model with bounded modeling errors is available. On the other side, the data-driven strategy directly computes the controller from the experimental data. Such techniques are also called direct methods and may be appealing in the cases such a control-oriented model is too time-consuming, too complex or too costly to obtain. The two strategies, model-based and data-driven, are complementary in sense that they do not address the same categories of problem.

Numerous direct methods have been proposed to try to achieve the best possible performance without using any plant model. Among them, some methods, like Iterative Feedback Tuning (IFT, [2]), Correlation-based-Tuning (CbT, [3]), Virtual Reference Feedback Tuning (VRFT, [4]) or Loewner Data-Driven Control (LDDC, [5]), can be designated as model-reference techniques. The principle of the model-reference problem is recalled on Figure 1. These approaches only require data from the plant $P$ and the desired closed-loop behaviour, given as a reference model transfer function $M\in\mathcal{RH}_{\infty}$. The objective is to design a controller that minimizes the error $\epsilon$ between the resulting closed-loop and the desired one $M$.

The choice of a reference model in data-driven control techniques is a critical step, see [1], [6], [7], [8]. Indeed, it should represent the desired closed-loop performances and be achievable by the plant at the same time. In [1], a method to build such a reference model, both reproducible by the system and having a desired behaviour, is proposed. It is applicable to Linear Time-Invariant (LTI) monovariable systems. The present paper aims at providing more applications of this method. The LDDC (Loewner Data Driven Control) algorithm is used to illustrate the impact of the choice of the reference model on the control design process.

This paper is an appendix of [1]. Its objective is to provide additional examples. This article is organized in five sections. Section 2 recalls the method proposed in [1]. Two control applications are given in Sections 3 and 4. The first considered plant is a continuous crystallizer. The second one is an open channel for hydroelectricty generation. Finally, conclusions and outlooks are proposed in Section 5.

2 Preliminary: Data-driven selection of a reference model

2.1 Definition of an achievable closed-loop

By achievable behaviour, we mean that the reference model corresponds to an internally stable closed-loop interconnection, see Figure 2. In other words, the ideal controller $K^{\star}$ should stabilize the plant internally. It is well known that the interconnection visible on Figure 2 is internally stable if and only if the sensitivity function $1-M$ is stable and if there are no compensation of instabilities in the open-loop, between the ideal controller $K^{\star}$ and the plant $P$.

Therefore, as explained in [1], the reference model $M$ is achievable by the plant if and only if $M$ is stable and satisfies:

[TABLE]

where $\left\{z_{i}\right\}_{i=1}^{n_{z}}$ and $\left\{p_{j}\right\}_{j=1}^{n_{p}}$ are respectively the unstable zeros and poles of the plant, which are assumed to be distinct. These interpolatory conditions are equivalent to having no compensation of instabilities between the ideal controller $K^{\star}$ and the plant $P$. When the plant has multiple RHP poles or zeros, derivative constraints must be added to (1), see [1] for further information.

2.2 Data-driven selection of an achievable reference model

The method proposed in [1] is briefly summed up in Algorithm 1. Details about the different steps are given in the following paragraphs.

2.2.1 Data-driven stability analysis

In order to define the interpolatory conditions of (1) for a given plant $P$, it is necessary to determine its RHP poles and zeros. This can be done in a data-driven way using the methods presented in [9] and [10].

The first step consists in a projection of the FRF measurements on the space $\mathcal{H}_{2}$ and $\overline{\mathcal{H}}_{2}$:

[TABLE]

where $P^{s}\in\mathcal{H}_{2}$ is its stable projection while $P^{as}\in\overline{\mathcal{H}}_{2}$ is its anti-stable part. On the basis of (2), it is possible to determine if the plant is stable or not. For further information on the projection method, see [10].

If the plant is unstable, it is then possible to estimate its RHP poles using its anti-stable projection $P^{as}$ as explained in [9]. These techniques are implemented in the PISA Toolbox [11].

2.2.2 Selection of an achievable reference model

As said earlier, an achievable reference model is stable and satisfies the interpolatory conditions given in (1), which are now well defined thanks to the estimation of the plant’s instabilities at the previous step.

The selection of a reference model depends on the nature of the plant. The proposed choice is detailed and justified in [1].

3 Application to a continuous crystallizer

The first considered application is the control of a continuous cooling crystallizer. This process is widely used in the chemical industry. It is a separation process which goal is to produce high-purity solids from liquids. The system is SISO: its input is the solute feed concentration $c_{f}(t)$ and its output is the solute concentration in the crystallizer $c(t)$. The state of the system is $x(t)=\left[n(L,t)\ c(t)\right]^{T}$, where $n(L,t)$ denotes the crystal size distribution. Physically, this system is described by population and mass balance equations. A complete mathematical model of this system is derived in [12].

The objective is to stabilize the plant around a desired steady-state $c(t)=c_{ss}=4.09mol/L$, which is just above the saturation concentration $c_{s}=4.038mol/L$, required for the crystals to be produced. For this steady state, as said in [13], the system is unstable and presents sustained oscillations which may degrade the crystals quality. Feedback control is therefore needed. This control problem has been treated in [13] through infinite-dimensional $\mathcal{H}_{\infty}$ synthesis, which is model-based, and in [14] thanks to a data-driven $\mathcal{H}_{\infty}$ synthesis.

When linearizing the system’s partial differential equations around the desired steady state, the crystallizer is characterized by an irrational transfer with an infinite number of poles, see [13], [12] or [14] for its expression. It is possible to evaluate numerically the frequency response of the linearized plant $P$ on a discrete frequency grid, see Figure 3. As in [14], a rational model $P_{502}$ of order 502 is obtained through a finite-difference method (see Figure 3). The poles and zeros of $P_{502}$ are given on Figure 4: the rational model is minimum phase and two unstable poles of value $3.83\times 10^{-5}\pm 0.848\times 10^{-2}\imath$ are visible.

In order to use the proposed method to select an achievable reference model, the considered frequency grid is much smaller: $N=500$ frequencies are considered, logspaced between $10^{-3}$ and 1 rad.s*-1*. The corresponding samples of the frequency response of the plant are estimated directly through the irrational transfer.

The first step of Algorithm 1 consists in the projection of the FRF measurements on the spaces $\mathcal{H}_{2}$ and $\overline{\mathcal{H}}_{2}$. The projection is given on Figure 5: the antistable projection fits the resonance while the stable part fits the rest of the frequency-response of the plant. Therefore the plant is unstable, as expected considering $P_{502}$ but also the previous studies of this system given in [13], [12] and [14].

As detailed in [9], the Hankel matrix of the antistable projection of the plant’s data is computed and a singular value decomposition is performed. The rank of this matrix gives us the order of the antistable projection, and consequently the number of RHP poles of the system $P$. The decomposition is visible on Figure 6: according to the drop after the second singular value, the system exhibits two unstable poles.

These two RHP poles are then estimated. Their value is given in Table 1. The obtained values are coherent with the ones found in [13] and with the RHP poles of the rational model $P_{502}$.

Remark 1

In [13], the RHP poles are estimated through a direct search method. This estimation is then used to factorize the plant’s expression to solve the mixed-sensitivity problem in the infinite-dimensional case.

Step 2 of Algorithm 1 consists in doing the same than in Step 1, but for the plant’s inverse, in order to determine whether the plant is minimum phase or not. The projection is then performed on the samples $\left\{\omega_{i},P(\imath\omega_{i})^{-1}\right\}$. The result is visible on Figure 7: the stable projection of the plant’s inverse fits the inverse of the plant’s frequency-response samples. Consequently, the plant is minimum phase.

Finally, an achievable reference model is selected according to Step 3 of Algorithm 1. The initial stable reference model $M$ is a first order transfer function:

[TABLE]

Since the plant is unstable and minimum-phase, the achievable reference model is chosen as $M_{f}=1-(1-M)B_{p}$, with $B_{p}$ defined according to the estimated RHP poles of the plant, see Table 1. The LDDC algorithm [15] is then applied. The identification of the controller is visible on Figure 8. The controller is reduced to a second-order model in order order to be compared with the one obtained in [14], its expression is given in (3).

[TABLE]

The finite-difference model is then used to simulate the closed-loop behaviour in time-domain. The identified controller is compared to the one obtained in [14], denoted $K_{[14]}$:

[TABLE]

The results are visible on Figure 11. The controller obtained in [14] allows to reach the desired steady-state faster with a less important overshoot. These better might be explained by two reasons: (i) model-reference control is really limitating when it comes to the expression of the closed-loop specifications and (ii) the reduction of the ideal controller degrades the closed-loop control performances.

First, let us investigate (i) the influence of the reduction of the ideal controller on the closed-loop control performances. As shown on Figure 9, the more the order of the identified controller is important, the more it will fit the frequency-response of the ideal controller. This is also visible on the time-domain simulations given on Figure 10, when simulating the passage to a new steady-state. A high order controller is more likely to give the desired closed-loop behaviour, specified by the reference model $M_{f}$.

Finally, let us investigate (ii) the influence of the specifications. In [14], the specifications are given as frequency weightings functions, giving more freedom to the desired closed-loop behaviour. To underline this aspect, the closed-loop $M_{14}$ reached by the controller $K_{[14]}$ obtained is [14] is taken as reference model. $M_{14}$ has been computed using the finite-difference model $P_{502}$. It leads to the identification of the second-order controller $K_{M_{14}}$:

[TABLE]

$K_{M_{14}}$ gives a better response time and a more important overshoot than $K_{[14]}$. However, the closed-loop performances induced by $K_{M_{14}}$ are much closer to the ones obtained by $K_{[14]}$ than the ones obtained using $M_{f}$ as a reference model.

To sum up, the first important aspect in order to obtain good closed-loop performances is the reduction of the ideal controller. The second one is that, for a given controller order, the choice of the reference model remains a critical step: choosing one desired transfer, even achievable, is limiting compared to robust specifications using frequency weightings.

To compare with a model-based method, in [13], an infinite-dimensional model-based $\mathcal{H}_{\infty}$ controller synthesis is performed on this very same application. The obtained controller is irrational and a reduction step is needed, which can be quite complicated. On this use-case, a reduced 8th order controller is obtained in [13]. On the other side, applying a structured model-based technique such as hinfstruct relies on an approximation of the irrational system, here the finite-difference model $P_{502}$, and is therefore time-consuming due to the complexity of the model.

For these reasons, data-driven control techniques are particularly indicated in this case. In [14], the considerations regarding the limitations due to the RHP poles and zeros of the plant are known through an initial stabilizing controller. Stability is guaranteed by the algorithm thanks to a test on the winding number. However, this method requires to build a fine frequency grid on which samples of the frequency-response of the plant are assumed to be available. Furthermore, the control design relies on iterative non-smooth optimization, which can be time-consuming and is sensitive to the considered initial stabilizing controller.

On the other side, the main strength of the LDDC algorithm is its simplicity. It is a one shot technique, it does not make strong assumptions and does not depend on an initial stabilizing controller.

4 Application to an hydroelectricity generation channel

The second application is an industrial problem provided by the French power producer EDF (Electricité de France). EDF uses water resources to generate green energy with run-of-the-river power plants. They rely on open-channel hydraulic systems that are non-linear and which dynamic depends on the operating point. Here, for simplicity we will consider one single operating point only.

Their physical model requires partial differential equations (namely Saint-Venant equations). In [16], a new irrational transfer function is proposed for open channels to represent the level-to-flow variations for any operating point. It is the solution of Saint-Venant equations under many assumptions. The system has two inputs, the entering and the outgoing flows $q_{e}$ and $q_{s}$, and one output, the water depth $h$. The transfer is given by:

[TABLE]

where

[TABLE]

where $x$ is the position of the measurement point on the channel, $Q_{0}$ the nominal flow, $L$ the length of the open channel. $B_{0}$, $\lambda_{1}(s)$ and $\lambda_{2}(s)$ depend on the canal configuration and the nominal flow (see [16]).

The system, which dynamic is visible in Figure 12, is extremely slow, has a delay behavior and a pole in limit of stability. Moreover, it has an infinite number of poles since the transfer function is irrational.

The system have been approximated in [16] by a 8th order rational transfer function depending on the nominal flow with input time delays: $\tau_{e}\simeq 500s$ and $\tau_{s}\simeq 1500s$, on $q_{e}$ and $q_{s}$ respectively. This approximation is shown in Figure 12 (dashed red). This model will be used to evaluate the performances of the identified controller.

The input flow $q_{e}$ in the open channel is seen as a disturbance (rain for example). The objective is to maintain the water depth to avoid flooding in the area. To this aim, the command signal is the output flow $q_{s}$. Therefore, only the transfer $G_{s}$ between the output flow $q_{s}$ and the water depth $h$ is considered, see (6).

In this example, the frequency approach is interesting since the system is represented by an irrational transfer function. Therefore, one cannot have a time-domain simulation. However we still can estimate samples of the frequency response of the system $\{\omega_{i},\Phi_{i}\},\ i=1\ldots N$, from which the ideal controller’s frequency response can be deduced. The samples of the frequency response $\Phi_{i}=G_{s}(\imath\omega_{i}),\ i=1\ldots N$ are extracted from the irrational transfer function $G_{s}$, for $N=500$ linearly spaced frequencies between $10^{-4}$ and $10^{-1}$ rad.s*-1*.

As for the previous application, the first step of the proposed method consist in determining the nature of the system thanks to the projection of its frequency-response measurements. The results of steps 1 and 2 of Algorithm 1 are visible on Figures 13 and 14 respectively.

Due to the presence of an integrator, the data of the plant is filtered by a bandpass filter for the stability analysis, see [10] for further explanations. Figure 13 shows that, except the integrator, the plant has no unstable poles. The mismatch between the projection and the plant’s data is due to the use of the bandpass filter. According to Figure 14, the plant is minimum phase: the stable projection fits the inverse of the plant’s data.

Therefore, the only constraint that the reference model should satisfy is $M(0)=1$, which would have been respected anyway to have zero tracking error. The objective is to stabilize the system and to obtain a faster dynamic. The reference closed-loop $M$ is chosen to be a second order continuous transfer function:

[TABLE]

with $\omega_{0}=10^{-4}rad.s^{-1}$ and $\xi=1$. It satisfies $M(0)=1$.

The frequency response of the ideal controller $K(\imath\omega_{i})$, which exactly provides the desired closed-loop behavior dictated by $M$ when placed in the closed-loop, is obtained as follows:

[TABLE]

The result of the identification step is given in Figure 15. The minimal realisation of the ideal controller is of order 137 and is stable. It is then reduced to a stable second order controller:

[TABLE]

Since we had no access to the EDF simulator, the 8th order rational transfer function of [16] is used to simulate the closed-loop with the 2nd order controller obtained by the LDDC. The results are shown in Figure 16: the resulting closed-loop achieve a response time of $4.84\times 10^{5}$s (134.4 hours) with no overshoot, while the system naturally has a response time of $2.47\times 10^{13}$s. The closed-loop dynamic is almost the objective one. The command signal is shown on Figure 17: it is reasonable, the maximum flow variation is around $8.3m^{3}.s^{-1}$, which is in the acceptable range for this application with a controller of order 2.

Finally, the performances in terms of disturbance rejection are shown on Figure 18: an input flow of $100m^{3}/s$ during four hours is considered. This disturbance would increase the water depth of $0.65m$ if not rejected. Figure 18 represents the tracking error when this disturbance is applied on the stabilized closed-loop: the water depth increases of $0.54m$ instead of the $0.65m$ without the controller. Finally, the disturbance is completely rejected $0.8\times 10^{5}s$ (2.3 hours) after its application. During the remaining time of the disturbance application, the disturbance does not affect the closed-loop system.

It should be noted that the closed-loop could be faster by increasing the frequency $\omega_{0}$ of the reference model $M$, see (7). However, taking $\omega_{0}=10^{-4}rad.s^{-1}$ (instead of $10^{-5}rad.s^{-1}$), leads to the identification of a second-order controller giving an oscillatory behaviour in closed-loop. The consequence is a significant overshoot, which is not acceptable in this application. Furthermore, the command signal $q_{s}$ rises to $80m^{3}/s$ when applying a step on the corresponding closed-loop.

In this application, the LDDC method is appealing since it does not require to simulate the complex system described by an irrational transfer function to obtain time-domain data. Only samples of the frequency response of the plant are needed, which can be estimated directly from the irrational transfer function. Moreover, one should notice that controlling such an infinite order model is also quite challenging even for model-based methods. A interesting perspective would be to try this controller on the EDF simulator instead of using the approximate model to validate the performances.

5 Conclusions

In this appendix, two additional examples illustrate the method proposed in [1] to select an achievable reference model for data-driven control purposes. The first one is a continuous crystallizer and the second one is an open-channel for hydroelectricty generation. These applications are representative of the class of systems for which data-driven control techniques is more appealing than model-based ones. Indeed, their model, which are irrational, are too complex for model-based control.

In both cases, the data-driven stability analysis proposed in [10] allows to draw the right conclusion concerning the nature of the plant. However, the second example considered in Section 4 shows that integrators are not handled by this technique. In Section 3, the continuous crystallizer is shown to be unstable and its unstable poles correspond to the ones previously found in the literature.

The techniques presented in [1] and recalled in Section 2 is then used to select an achievable reference model. The resulting ideal controllers are stable and, in the case of the continuous crystallizer, do not compensate the plant’s instabilities. The LDDC algorithm is then used to identify reduced-order controllers.

As shown by the application of the LDDC algorithm on the continuous crystallizer, the proposed method does not compensate the fact that choosing a reference model constitutes a very limited specifications requirement compared to robust specifications. Further work will treat the choice of a reference model in the multivariable case: the problem becomes a tangential interpolation one. Another outlook would be to investigate the possibility to tune the reference model to have a better control on the corresponding closed-loop performances.