Data-based stabilization of unknown bilinear systems with guaranteed basin of attraction

Motivated by the goal of having a building block in the direct design of data-driven controllers for nonlinear systems, we show how, for an unknown discrete-time bilinear system, the data collected in an offline open-loop experiment enable us to design a feedback controller and provide a guaranteed under-approximation of its basin of attraction. Both can be obtained by solving a linear matrix inequality for a fixed scalar parameter, and possibly iterating on different values of that parameter. The results of this data-based approach are compared with the ideal case when the model is known perfectly.


Introduction
Direct data-driven control design aims at obtaining a control law based only on input-output data collected from the system during an experiment, thereby avoiding altogether the identification of a model of the system from the data. Within the literature of direct data-driven control, notable approaches are iterative feedback tuning [20], virtual-reference feedback tuning [9], iterative correlationbased tuning [23] to name a few. When data are assumed to be generated by an underlying linear system, a number of approaches are available [33,18,2,30,13,38]. However, sensibly fewer address the case when data are generated by an underlying nonlinear system, see [17,10,28].
To address the intrinsic difficulty of dealing with the control design of unknown nonlinear systems, a natural approach is to reduce their complexity by considering the system evolution along a given Lyapunov function. This classical control theoretic analysis is enhanced by nonparametric regression methods from machine learning to cope with the large uncertainty in the model [4] and is performed using a sufficiently dense set of samples taken from the system. Analytical guarantees of stability and safety are then obtained relying on additional tools from robust control and optimization [41]. The approach of [35,17] to reduce the complexity of controlling unknown nonlinear systems consists of considering systems with a well-defined relative degree, in such a way that the uncertainty only appears in the form of two Lie derivatives of the output function along the system vector fields. Once the dynamics has been discretized, the key observation from sampleddata control theory is that these uncertain functions are constant between sampling times for a sufficiently high sampling rate.
A different approach to data-driven control of nonlinear systems has been recently taken in a series of works that use the nonparametric representation of dynamical systems via Hankel matrices of finite-size input-output data proposed in [42]. On one hand, this representation has given rise to data-enabled predictive controllers where the effect of the nonlinearity is taken into account by a regularized optimization problem [12,22]. On the other hand, it inspired a data-dependent parametrization of the closedloop system that reduces the control design to semidefinite programs where the nonlinearity is dealt with as a process disturbance [15]. Further results along this research thread have been proposed in [3]. While these results make possible to deal with nonlinear systems, they provide local stability results. Very recently, within the research thread of [42,Thm. 1], there have been efforts to go beyond the local nature of the results for special classes of nonlinear systems, studying data-driven control of second-order discrete Volterra systems [31] and polynomial systems [19].
The goal of this paper is to characterize another notable class of nonlinear systems for which nonlocal datadriven control results can be established, namely bilinear systems. The reason for focusing on bilinear systems is threefold. In spite of their simple nonlinear structure, applying Carleman linearization to a generic continuous-time input-affine nonlinear system yields a continuous-time bilinear system with a larger state plus a remainder (see [7,25]), so bilinear systems can be used as universal approximators of input-affine nonlinear systems [32, p. 110]. This last consideration specifically motivates the proposed arXiv:2004.11630v1 [eess.SY] 24 Apr 2020 data-driven control scheme for bilinear systems, which is envisioned to be a building block in future work on direct data-driven control of input-affine nonlinear systems (see also the discussion in Remark 3). A second motivation is to provide a method alternative to sum-of-squares programming for polynomial control systems [19] to directly design data-driven controllers of bilinear systems. Finally, bilinear systems are interesting per se as meaningful models for a number of relevant applications in engineering, biology and ecology [27,8].
Many model-based approaches have been proposed for control of bilinear systems such as [5,37,1,24], and we refer the reader to [24, §1] for a thorough overview. Such model-based approaches assume the knowledge of the parameters of the bilinear system. When these are not known from first-principles considerations, one can resort to system identification techniques tailored for bilinear systems, and then apply one of the model-based approaches above. Some of these indirect data-driven methods for system identification are [16,11,34], see also [40,Part II] for an overview. Although combining the aforementioned system identification techniques with model-based design constitutes a natural and valid way to control a bilinear system, we aim here at exploring the less-investigated direct control design of a bilinear system based on data (avoiding altogether a system identification step generally nontrivial in a nonlinear setting). We show that under mild assumptions (see Assumption 1 below), it is indeed possible to design stabilizing control policies directly from data. We also show via simulations that our approach compares well with a model-based design that has perfect knowledge of the parameters of the system, regardless of whether this knowledge derives from first-principles considerations or from a preliminary system identification step.
In the case of data generated by an underlying linear system, the fundamental result [42,Thm. 1] has been shown in [15] to allow direct data-driven design of (optimal) feedback controllers (with robustness to noise) for linear systems through linear matrix inequalities (LMI) [6] and the local stabilization of nonlinear systems through semidefinite programs. In the case of data generated by an underlying bilinear system, the arguments in [15] need substantial modifications to counteract the nonlinear term appearing in the bilinear system and to explicitly provide an estimate of the region of attraction. Thus, we need to resort to tools from robust control (such as [29,24], see Fact 1 below) besides more standard ones from linear matrix inequalities. Some conservatism is introduced in these steps compared to a model-based approach, as illustrated in Section 4. Similarly to the model-based approaches [24,37] and, partially, to [5,1], we also adopt a linear state feedback and a quadratic Lyapunov function in the design of the closed-loop system. Alternatives are based on rational polynomial controllers and sum-of-squares programming [39]. The choice of linear controllers is restrictive compared to nonlinear state feedback (and the actual basin of attraction has not an ellipsoidal shape), but are dictated by the the desire of obtaining a computationally tractable result in the form of linear matrix inequalities (after fixing a scalar parameter). However, the main difference with those model-based approaches is that we design here the linear state feedback and the quadratic Lyapunov function without relying on the knowledge of the bilinear system matrices, which we aim to substitute instead through data collected from the bilinear system.
Tuning a feedback controller based only on a limited number of open-loop data, which gives a guaranteed subset of the basin of attraction for a bilinear system, is the main contribution of this paper.
Structure. The considered problem is formulated in Section 2. In Section 3 we provide our data-based controller for the unknown bilinear system with a guaranteed under-approximation of its basin of attraction, as a main result. Section 4 compares this solution with a modelbased one on a numerical example.
Notation. For a matrix A, A denotes the induced 2-norm. For a symmetric matrix A B B C , we may use the shorthand writing [ A B C ]. I denotes an identity matrix of appropriate dimensions.

System description and problem formulation
Consider the discrete-time bilinear system where x ∈ R n is the state, u ∈ R is the input, and the system matrices have dimensions A ∈ R n×n , B ∈ R n , D ∈ R n×n . Our choice to consider a scalar input in (1) is motivated in Remark 2 after we have outlined our approach. The matrices A, B, D are completely unknown apart from a bound on the matrix norm of D as follows.
Our objective is to design a controller u = Kx for the bilinear system in (1) based only on data collected from an off-line experiment (namely, without identifying the matrices A, B, D) and give a guaranteed under-approximation of the basin of attraction of the origin for the closed-loop system. The off-line experiment of duration T (with T > 0) collects the input and state sequences u(0), u(1), . . . , u(T − 1) and x(0), x(1), . . . , x(T ). These are organized as and allow computing the auxiliary quantity Following [15], we reparametrize the gain K by a matrix G K and give in the next lemma an equivalent representation of (1) in closed loop with u = Kx, which depends on data, except for the matrix D.
Then, system (1) with state feedback u = Kx and K = U 0,T G K has the equivalent representation This closed-loop matrix is, by (3), since the data in (2) satisfy X 1,T = AX 0,T + BU 0,T + DV 0,T , and this proves the statement.
The reparametrization G K is a decision variable that we tune to achieve our control objective. Based on G K and on data, we define for compactness so that the closed-loop representation in (4) becomes where D is highlighted. As mentioned earlier, we aim at giving a guaranteed under-approximation of the basin of attraction of the closed-loop system in (6). We do so by considering a quadratic Lyapunov function with Q = Q 0 and imposing the strict decrease of V g D (x)x − V (x) for the dynamics in (6). The last quantity is easily computed as in the next claim.
Note that for D = 0, (1) becomes linear and (8) reduces to N D (x) = A c QA c − Q, corresponding to the classical Lyapunov condition for discrete-time linear systems. We by designing the decision variables G K , which determines K = U 0,T G K , and Q, which will be optimized so that the volume of the ellipsoid E Q is maximized. The design will be based only on data, and return the ellipsoid E Q as a guaranteed under-approximation of the basin of attraction. We summarize the system description and control objective illustrated in this section as follows.
Problem 1. Based only on data collected from an off-line experiment as in the matrices in (2), obtain a controller u = Kx for (1) such that for the closed-loop system, the origin has a guaranteed basin of attraction. The data-based design is performed based on the decision variables G K and Q of the quadratic Lyapunov function in (7), a sublevel set of which gives the guaranteed basin of attraction.
Some remarks are in order.
is related to the "quality" of the experimental data. In fact, condition (3) expresses the property that the data are sufficiently rich so that the system dynamics can be parametrized directly in terms of the matrices in (2). A key property established in [42] is that, for linear systems, X 0,T is full-row rank (thus, a solution G K to (3) exists) when the experiment is carried out using a sufficiently exciting input signal. An extension of this property to nonlinear systems is discussed in [14] where it is shown that under prior knowledge of an upper bound on the nonlinearity (in fact, on D in the present case of bilinear systems) one can always design experiments so that (3) is fulfilled.
Remark 2. The present analysis can be extended to bilinear systems with input u ∈ R m and m ≥ 2. For m = 2, (1) can be written for u = (u 1 , u 2 ) as We can define U This expression shows by comparison with (4) that the case for m = 2 can be treated using the same procedure we develop in the presence of a single unknown D, and this consideration easily generalizes to m larger than 2. For this reason we focus on the essential case with input u ∈ R. Remark 3. The universal approximation property of bilinear systems mentioned in Section 1 holds with respect to continuous-time nonlinear systems. We focus here on discrete-time bilinear systems since the data in (2) are samples obtained from experiments. However, the same results would apply to continuous-time bilinear systems if X 1,T in (2c) is replaced by the samples of the state derivatives. These results would then lend themselves to the analysis of a bilinear approximation of continuous-time nonlinear systems (provided disturbances are accounted for).

Data-based solution with guaranteed basin of attraction
In Section 2, we showed that data allow expressing (1) in closed loop with u = Kx as (6) (by introducing the reparametrization G K of K). Data, however, did not allow us to completely remove the matrices of model (1). In particular, g D (x) in (6) still contains two instances of matrix D (namely, DxK and FDH), which can both be interpreted as perturbation of the matrix A c . In this section we first address the former, which is more standard and occurs analogously for model-based design of a bilinear system (see, e.g., [24]), and then the latter, which is motivated by our desire to solve Problem 1 based only on data and calls for the matrix norm bound in Assumption 1.
Before presenting the developments of this section, we recall an auxiliary result from [29], which has been reported in a convenient form as [24,Lemma 1] and is related to the S-procedure [6, §2.6.3]. In particular, [24,Lemma 1] implies the next fact.
if there exists a scalar e such that Remark 4.
[1] considers a model-based setting similar to [24] returning a polytope as a subset of the basin of attraction. However, since [1, Thm. 2] does so by dilating the polytope and including in that a quadratic sublevel set, we refer to the naturally quadratic approach in [24].
With Fact 1 we are in a position to develop this section. The next lemma addresses the term DxK in g D (x) in (6). Specifically, it shows that as long as we restrict the analysis to a sublevel set E Q (defined in (9)) of the Lyapunov function V in (7) (where Q itself is a decision variable determining the size of this sublevel set), strict decrease of V along solutions is guaranteed (N D (x) ≺ 0) since N D (x) determines V g D (x)x − V (x) as in Lemma 2.
Proof. The proof follows closely [24], but is reported for self-containedness. Define for compactness and note for the following that (13) implies Q 0 and τ > 0. By Schur's complement (with respect to lowest block −Q −1 ), (13) is equivalent, by (14), to By Schur's complement, this inequality is equivalent to Rearranging rows and columns of this inequality gives which is equivalent to (13). We want to put this inequality in a form where we can apply Fact 1. Then, we pre-and postmultiply the previous inequality by the block diagonal matrix with entries I, I, (Q 1/2 ) −1 , I (where Q 1/2 is the unique symmetric, positive definite square root matrix for Q = Q 0 [21, Thm. 7.2.6], so that Q = Q 1/2 Q 1/2 ) and apply Schur's complement (with respect to the lowest block − 1 τ I) to obtain with some computations  Note that x Qx = (x Q 1/2 )(Q 1/2 x), hence for all x such that x Qx ≤ 1, x Q 1/2 ≤ 1. With this observation and by Fact 1 we conclude (after some simplifications) that for all x such that x Qx ≤ 1. We show now that this is equivalent to the conclusion of the lemma. Define for compactness so that (15) is equivalent, after some computations, to By Schur's complement, we obtain that which is equivalent by (8) The next lemma addresses the term FDH in g D (x) in (6). Specifically, it shows that as long as the matrix D is bounded in norm by δ as in Assumption 1, we can obtain a matrix inequality depending only on δ and guarantee that Lemma 3 and its conclusions hold for all such D, which is key to obtain a fully data-based solution to our problem.  (16) then (13) holds.
Proof. Note that from F = I in (5), (13) is equivalent to and this equation has the same structure as G + MDN + ND M ≺ 0 in Fact 1, since D ≤ δ (δ > 0) by Assumption 1. Indeed, by making the suitable correspondences between the quantities of this lemma and those of Fact 1, the existence of 2 such that (16) holds (corresponding to (12) of Fact 1) guarantees that (13) (corresponding to (11) of Fact 1) holds for D as in Assumption 1.
Lemma 4 enables us to generalize the conclusions of Lemma 3 for all D with D ≤ δ, so that we do not need to rely on the knowledge of D (as it would be the case in a model-based scheme), but just on its (possibly loose) norm bound δ. The matrix inequality (16) of Lemma 4 (where only δ appears), however, contains products of decision variables and inverses of decision variables. We address this in the next proposition, which obtains a matrix inequality that is as close as possible to an LMI (hence efficient to solve) and expresses explicitly the matrix inequality in terms of the available data. This proposition is the main result of this paper. Proposition 1. Under Assumption 1, suppose there exist 1 ∈ R, 2 ∈ R, Y ∈ R n×T and P = P ∈ R n×n such that and set Q = P −1 , G K = Y P −1 . Then, (i) for the dynamics in (6) corresponding to D, the Lya- (ii) the origin is asymptotically stable for (1) with con- x and its basin of attraction contains the set E Q .
Proof. We begin showing that inequalities (16) and (17a) are equivalent, noting for the following that (17a) implies P 0. With the definitions in (5), (16) is equivalent to By pre-and postmultipling this inequality by the block diagonal matrix with entries Q −1 , Q −1 , I, I, I and by setting Q = P −1 , G K = Y P −1 as in the statement of the proposition, the last inequality is equivalent to To avoid the simultaneous presence of τ and 1/τ , this inequality is equivalent to the next one by pre-and postmultiplying by the block diagonal matrix with entries I, 1 τ I, I, I, I and setting 1 = 1/τ : which is exactly (17a). After these manipulations, the conclusions of the proposition follow readily. Indeed, the fact that (17a) holds, implies that (16) holds, and then by Lemmas 4 and 3 that D as in Assumption 1 satisfies N D (x) ≺ 0 for all x ∈ E Q . By Lemma 2, (i) follows. (17b), which is equivalent to I = X 0,T G K , and Lemma 1 ensure that (4) or, equivalently, (6) are an equivalent representation of (1) with controller u = Kx = U 0,T G K x. Standard Lyapunov theorems give then (ii). Proposition 1 effectively solves Problem 1. Indeed, if a solution to (17) is found (which is based on data from an off-line experiment), then we have a controller K and a guaranteed basin of attraction in terms of the set E Q .
The matrix inequality (17a) in Proposition 1 is convenient because, after fixing the scalar 1 , it is an LMI in the decision variables 2 , Y , P . A line search with respect to 1 on top of solving this LMI is typically preferable than solving directly the bilinear matrix inequality in (17a). Note also that model-based approaches for controlling bilinear systems encounter such a situation, and fix one of the parameters directly [24] or in an iterative way [37].
A conclusion of Proposition 1 is that the basin of attraction of the origin contains the set E Q = E P −1 . It is quite natural to maximize the volume of this ellipsoid, which is proportional to det(P ), as is done in the modelbased setting of [24]. (Other size criteria can be optimized, see the discussion in [36, §2.2.5.1].) This leads to the next immediate corollary. Corollary 1. Let Assumption 1 hold. If there exist a solution to the next optimization problem in the decision variables 1 ∈ R, 2 ∈ R, Y ∈ R n×T and P = P ∈ R n×n minimize − log det(P ) subject to (17a), (17b), then the conclusion of Proposition 1 holds.
Finally, since we are considering a quadratic Lyapunov function and as is done in the model-based solutions [24,37], the very same arguments leading to Proposition 1 yield exponential (instead of asymptotic) stability by strengthening a little the matrix inequality in (17a). This is stated in the next corollary, whose proof is thus omitted.
Corollary 2. For µ ∈ (0, 1), suppose that the assumptions of Proposition 1 can be satisfied after replacing the element (1, 1) of the matrix in (17a) (i.e., −P ) with −µP . Then, (i) for the dynamics in (6) corresponding to D, the Lya- (ii) the origin is exponentially stable for (1) with controller u = Kx = U 0,T G K x = U 0,T Y (X 0,T Y ) −1 x and its basin of attraction contains the set E Q .

Numerical example
We consider for (1) the matrices  which are taken from [5, §5]. Our design does not rely on their knowledge, but simply on the data generated according to them and a bound δ of D . In particular, we consider δ = 0.7637, which overapproximates by 20% the actual D = 0.6364. Moreover, we will use the matrices in (18) to compare our data-based design with a modelbased design in Section 4.2. We note that the comparison is made with a model-based design that has perfect knowledge of the parameters of the system. Getting to perfectly know the parameters would correspond to the ideal case even for a preliminary system identification step. We show in this section that our designed controller performs comparably to such a model-based design, in spite of being tuned only on an offline experiment. We consider T = 10. In Figure 1, we show the input and state sequences giving (2) and generated according to the matrices in (18). We note that A being unstable is challenging because a suitable control action has to be designed to modify by feedback the system evolution in a neighborhood of the origin (without the "help" of a stable linear part) and the generated data quickly grow large as shown in Figure 1, thereby impacting the numerical accuracy of the procedure.

Data-based solution
By using Corollary 1, the data-based solution implemented in this section is as follows. In particular, we opt for fixing the scalar variable 1 , solve an LMI, and perform a line search on 1 .
2. We solve the next optimization problem in the decision variables 2 ∈ R, Y ∈ R n×T and P = P ∈ R n×n minimize − log det(P ) which is an LMI. By denoting the solution P =: P DB , we then obtain G K = Y P −1 DB and the controller gain as K DB := U 0,T G K .
3. We iterate on the selection of 1 in case of, e.g., infeasibility.
We implement this scheme (and the model-based one in Section 4.2) through YALMIP [26]. For a value of 1 = 0.8, we obtain The evolution of x when u = K DB x is used in (1) is given in Figure 2 in the top plot as a phase portrait (solid colored lines) and in the middle plot as a time evolution.

Comparison with model-based solution
For (1) with matrices in (18), we use the model-based solution in [24] for comparison. This model-based solution is also not an LMI, unless the scalar parameter 1 is fixed (as in the data-based solution) and a line search is performed.
2. We solve the optimization problem in the decision variables y ∈ R n and P = P ∈ R n×n minimize − log det(P ) which is an LMI. By denoting the solution P =: P MB , we then obtain the controller gain as K MB := y P −1 MB .
3. We iterate on the selection of 1 in case of, e.g., infeasibility.
For the same value of 1 = 0.8 as in Section 4.1, we obtain The evolution of x when u = K MB x is used in (1) is given in Figure 2 in the top plot as a phase portrait (dotted colored lines) and in the bottom plot as a time evolution.

Discussion
Finally, we compare the performance of the data-based solution against the model-based solution by performing a thorough line search on the parameter 1 , which we fixed before in order to be able to solve an LMI. The result is in Figure 3. Only values of 1 where an optimal solution was returned by YALMIP, are displayed (in particular, this did not happen for the model-based solution with values of 1 between 0.2 and 0.4).
The top plot represents the determinants of the matrices P DB and P MB , which was considered since it is proportional to the volume of the ellipsoids that are guaranteed to be in the basin of attraction of the closed-loop system.
In the middle plot, the logarithms of these determinants are also provided since they are the actual objective functions in the optimization problems of Sections 4.1-4.2.
As expected, the model-based solution provides ellipsoids with larger sizes (e.g., det(P MB ) = 60.03 for 1 = 0.4). For the given example, it appears from Figure 3 that the data-based solution performs better for small 1 , whereas it performs worse than the model-based solution for large 1 . We note that log det is actually more representative of the actual difference between the two solutions. Indeed, for values of 1 around 1, the two solutions are not so distant, as is confirmed by the illustration of Figure 2 where the corresponding ellipsoids are also depicted in the top plot (solid and dotted black curves).
Finally, since the model-based and data-based solution share a similar structure, we expect that they lead to similar feedback gains. This is indeed confirmed by their relative difference in norm in the bottom plot in Figure 3.
In summary, our designed controller presents in these simulations a similar performance to the model-based design, where the former relies on an offline experiment and the latter on the perfect knowledge of system parameters.

Conclusions
We proposed a direct data-driven design for bilinear systems, which comes with a guaranteed subset of the basin of attraction.
The main goal of future work is applying this scheme as a building block for data-driven control of input-affine nonlinear systems (by approximating the latter through Carleman linearization). Closely related topics of future work are a study of the modifications needed to cope with noisy data, and of the tradeoffs with schemes based on sum-of-squares programming for bilinear systems.