Globally stabilizing state feedback control design for Lotka-Volterra systems based on underlying linear dynamics∗

Preprints, 1st IFAC Conference on Modelling, Identification and Preprints, IFAC Conference Modelling, Identification and Control of 1st Nonlinear Systems on Preprints, 1st IFAC Conference Conference on Modelling, Modelling, Identification Identification and and Preprints, IFAC on Control of 1st Nonlinear Systems Available online at www.sciencedirect.com June 24-26, 2015. Saint Petersburg, Russia Control of Nonlinear Systems Control of Nonlinear Systems June 24-26, 2015. Saint Petersburg, Russia June June 24-26, 24-26, 2015. 2015. Saint Saint Petersburg, Petersburg, Russia Russia

ScienceDirect

IFAC-PapersOnLine 48-11 (2015) 1000–1005

Globally stabilizing state feedback control Globally stabilizing state feedback control Globally stabilizing state feedback control design for Lotka-Volterra systems based on design for Lotka-Volterra systems based on design for Lotka-Volterra systems based on  underlying linear dynamics   underlying underlying linear linear dynamics dynamics Attila Magyar ∗∗ Katalin M. Hangos ∗,∗∗ ∗,∗∗ Attila ∗,∗∗ Attila Magyar Magyar ∗∗ Katalin Katalin M. M. Hangos Hangos ∗,∗∗ Attila Magyar Katalin M. Hangos ∗ ∗ Department of Electrical Engineering and Information Systems, ∗ Department of Electrical Engineering and Information Systems, ∗ Department of Electrical Electrical Engineering and Information Systems, University of Pannonia Egyetem street 10.and Veszpr´ em, H-8200 Hungary Department of Engineering Information Systems, University of Pannonia Egyetem street 10. Veszpr´ eem, H-8200 Hungary University Pannonia Egyetem street 10. Veszpr´ m, H-8200 Hungary (e-mail:of {magyar.attila; hangos.katalin}@virt.uni-pannon.hu). University of Pannonia Egyetem street 10. Veszpr´ e m, H-8200 Hungary {magyar.attila; hangos.katalin}@virt.uni-pannon.hu). ∗∗ (e-mail: (e-mail: {magyar.attila; hangos.katalin}@virt.uni-pannon.hu). Process Control Research Group, Control and Automation Research ∗∗ (e-mail: {magyar.attila; hangos.katalin}@virt.uni-pannon.hu). ∗∗ Process Control Research Group, Control and Automation Research ∗∗ Process Control Research Group, Control and Automation Automation Research Institute, HAS, Research Kende street 13-17. Budapest, H-1111, Hungary Process Control Group, Control and Research Institute, HAS, Kende street 13-17. Budapest, H-1111, Hungary Institute, Institute, HAS, HAS, Kende Kende street street 13-17. 13-17. Budapest, Budapest, H-1111, H-1111, Hungary Hungary Abstract: Using translated X-factorable phase space transformations and nonlinear variable Abstract: Using translated X-factorable phase space transformations and nonlinear variable Abstract: Using X-factorable phase space transformations nonlinear variable transformations dynamically similar linear ODE model is associated and to the Lotka-Volterra Abstract: Usingaa translated translated X-factorable phase space transformations and nonlinear variable transformations dynamically similar linear ODE model is associated to the Lotka-Volterra transformations a dynamically similar linear ODE model is associated to the Lotka-Volterra system models with a positive equilibrium point. This enables to use a linear full state feedback transformations a dynamically similar linear ODE model is associated to the Lotka-Volterra system with equilibrium point. enables to aa linear full state feedback system models models with a a positive positive equilibrium point. This This enables to use use the linear full state feedback controller to stabilize the system in controllable cases, that leaves openfull loop equilibrium system models with a positive equilibrium point. This enables to use a linear state feedback controller to stabilize the system in controllable cases, that leaves the open loop equilibrium controller to stabilize stabilize the system system in controllable controllable cases, thatproblem leaves the the open loopcase equilibrium point unchanged. The linear state feedback controller design in the general has also controller to the in cases, that leaves open loop equilibrium point unchanged. The linear state feedback controller design problem in the general case has also point unchanged. The linear state feedback controller design problem in the general also been formulated to ensure the compartmental property of the closed loop systemcase fromhas which point unchanged. The linear state feedback controller design problem in the general case has also been formulated to ensure the compartmental property of the closed loop system from which been formulated to ensure the compartmental property of the closed loop system from which the existence of a diagonal Lyapunov function follows. Further extension has been obtained by been formulated to ensure the compartmental property of the closed loop system from which the of a function Further extension obtained the existence existence ofre-parametrization a diagonal diagonal Lyapunov Lyapunov function follows. follows. Further extension has has been been obtained by by using the timeof transformation defined for quasi-polynomial models. the existence a diagonal Lyapunov function follows. Further extension has been obtained by using the time re-parametrization transformation defined for quasi-polynomial models. using using the the time time re-parametrization re-parametrization transformation transformation defined defined for for quasi-polynomial quasi-polynomial models. models. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: quasi-polynomial systems, Lotka-Volterra models Keywords: quasi-polynomial quasi-polynomial systems, Lotka-Volterra Lotka-Volterra models Keywords: Keywords: quasi-polynomial systems, systems, Lotka-Volterra models models 1. INTRODUCTION power of such systems. Recently, the positive polynomial 1. INTRODUCTION INTRODUCTION power such Recently, the positive polynomial 1. power of ofclass suchofsystems. systems. Recently, the positive polynomial system chemicalRecently, reaction the networks (CRNs) have 1. INTRODUCTION power of such systems. positive polynomial system class of chemical reaction networks (CRNs) have system class of to chemical reaction networks networks (CRNs) have A wide range of systems can only be tackled using nonlin- system also been used design stabilizing polynomial feedbacks class of chemical reaction (CRNs) have A wide wide range range of of systems systems can can only only be be tackled tackled using using nonlinnonlin- also been used to design stabilizing polynomial feedbacks A also been used to design stabilizing polynomial feedbacks ear techniques (Isidori (1995)). Most of themusing are applicain Szederk´ enyi et al. (2013). A wide range of systems can only be tackled nonlinalso been used to design stabilizing polynomial feedbacks ear techniques techniques (Isidori (Isidori (1995)). (1995)). Most Most of of them them are are applicaapplica- in Szederk´ eenyi et al. (2013). ear Szederk´ al. (2013). ble only for a narrow of nonlinear the in ear techniques (Isidoriclass (1995)). Most of systems, them arewhile applicain Szederk´ enyi nyi et et al. (2013). ble only for a narrow class of nonlinear systems, while the Previous work in the field of quasi-polynomial systems ble only for aa narrow class of nonlinear systems, while the more generally applicable methods suffer from computaPrevious work in the field of quasi-polynomial systems ble only for narrow class of nonlinear systems, while the Previous work in the field of more generally applicable methods suffer from computa(Figueiredo et al. (2000)) that the globalsystems stabilPrevious work in the fieldproved of quasi-polynomial quasi-polynomial systems more generally applicable suffer from computational complexity problems.methods One possible way of balancing (Figueiredo et al. (2000)) proved that the global stabilmore generally applicable methods suffer from computa(Figueiredo et al. (2000)) proved that the global stabiltional complexity problems. One possible way of balancing ity analysis of quasi-polynomial systems is equivalent to (Figueiredo et al. (2000)) proved that the global stabiltional complexity problems. One possible way of balancing between general applicability andpossible computational feasibil- ity analysis of quasi-polynomial systems is equivalent to tional complexity problems. One way of balancing ity analysis of quasi-polynomial systems is equivalent to between general applicability and computational feasibilthe feasibility of a linear matrix inequality (LMI). This ity analysis of quasi-polynomial systems is equivalent to between general applicability and computational feasibility is to find nonlinear system classes with good descriptive the feasibility of a linear matrix inequality (LMI). This between general applicability and computational feasibilthe feasibility of a linear matrix inequality (LMI). This ity is to find nonlinear system classes with good descriptive has lead to the conclusion in Magyar et al. (2008), that the feasibility of a linear matrix inequality (LMI). This ity is find nonlinear system with descriptive power characterized structure, and utilize this has lead to the conclusion in Magyar et al. (2008), that ity is to tobut findwell nonlinear system classes classes with good good descriptive lead conclusion in et power but well characterized structure, and utilize utilize this has the globally stabilizing nonlinear state feedback designthat for has lead to to the the conclusion in Magyar Magyar et al. al. (2008), (2008), that power but well characterized structure, and this structure when developing control design methods. the globally stabilizing nonlinear state feedback design for power but well characterized structure, and utilize this the globally stabilizing nonlinear state feedback design for structure when when developing developing control control design design methods. methods. quasi-polynomial systems is equivalent to a bilinear matrix the globally stabilizing nonlinear state feedback design for structure quasi-polynomial systems is equivalent to a bilinear matrix structure when developing control design methods. quasi-polynomial systems is equivalent equivalent to a bilinear bilinear matrix The class of quasi-polynomial (QP) systems plays an quasi-polynomial inequality. Although the solution of a to bilinear matrix insystems is a matrix The class class of of quasi-polynomial quasi-polynomial (QP) (QP) systems systems plays plays an an inequality. Although the solution of matrix inThe inequality. Although the solution of aaa bilinear bilinear matrix inimportant role in the theory of nonlinear and nonnegative equality is an NP hardthe problem an iterative LMImatrix algorithm The class of quasi-polynomial (QP) systems plays an inequality. Although solution of bilinear inimportant role role in in the the theory theory of of nonlinear nonlinear and and nonnegative nonnegative equality is an NP hard problem an iterative LMI algorithm important equality is an NP hard problem an iterative LMI algorithm dynamical systems because nonlinear systems with smooth equality could be isused (VanAntwerp andanBraatz (2000)). Recently, important role in the theory of nonlinear and nonnegative an NP hard problem iterative LMI algorithm dynamical systems systems because because nonlinear nonlinear systems systems with with smooth smooth could be used (VanAntwerp and Braatz (2000)). Recently, dynamical could be (VanAntwerp and Braatz Recently, nonlinearities can because be transformed quasi-polynomial methods have been for stabilizing dynamical systems nonlinearinto systems with smooth improved could be used used (VanAntwerp and developed Braatz (2000)). (2000)). Recently, nonlinearities can be transformed into quasi-polynomial improved methods have been developed for nonlinearities can be transformed quasi-polynomial improveddesign methods have been based developed for stabilizing stabilizing form (Hern´ andez-Bermejo and Fair´eninto (1995)). This means, improved feedback of QP systems on control Lyapunov nonlinearities can be transformed into quasi-polynomial methods have been developed for stabilizing form (Hern´ a ndez-Bermejo and Fair´ e n (1995)). This means, feedback design of QP systems based on control Lyapunov form (Hern´ aandez-Bermejo and een means, feedback design design of QP QP and systems based on control control Lyapunov that any applicable method forFair´ quasi-polynomial functions (see Magyar Hangos (2013)), and Lyapunov also using form (Hern´ ndez-Bermejo and n (1995)). (1995)). This Thissystems means, feedback of systems based on that any any applicable method forFair´ quasi-polynomial systems functions (see Magyar and Hangos (2013)), and also using that applicable for quasi-polynomial systems functions (see Magyar Hangos (2013)), also can be regarded asmethod a general technique for nonlinear the underlying reducedand linear dynamics in and Magyar et al. that any applicable method for quasi-polynomial systems functions (see Magyar and Hangos (2013)), and also using using can be regarded as a general technique for nonlinear the underlying reduced linear dynamics in Magyar et al. can be regarded as aa general technique for nonlinear the underlying reduced linear dynamics in Magyar al. systems. (2013). All of the above attempts, however, have et used can be regarded as general technique for nonlinear the underlying reduced linear dynamics in Magyar et al. systems. (2013). All of the above attempts, however, have used systems. (2013). All All nonlinear of the the above above attempts, however, have used polynomial feedback to achieve their control systems. (2013). of attempts, however, have used polynomial nonlinear feedback to achieve their control QP-systems are invariant under quasi-monomial transfor- polynomial nonlinear feedback to achieve their control goals. QP-systems are are invariant invariant under under quasi-monomial quasi-monomial transfortransfor- polynomial nonlinear feedback to achieve their control QP-systems mation (Figueiredo et al. (2000)), this enables to partition QP-systems are invariant under quasi-monomial transfor- goals. goals. mation (Figueiredo et al. (2000)), this enables to partition goals. mation (Figueiredo et (2000)), this to them into equivalence Lotka- The aim of this paper is to apply a linear state feedback mation (Figueiredo et al. al.classes (2000)),represented this enables enablesby to aapartition partition them into into equivalence classes represented by Lotka- The aim of this paper is to apply aa linear state feedback them equivalence classes represented aa Lotkaaim paper is apply feedback Volterra (LV) system. Dynamic similarity canby also be de- The structure forthis a special Lotka-Volterra systems usthem into equivalence classes represented by LotkaThe aim of of this papersubset is to to of apply a linear linear state state feedback Volterra (LV) system. Dynamic similarity can also be destructure for a special subset of Lotka-Volterra systems usVolterra (LV) system. Dynamic similarity can also be destructure for a special subset of Lotka-Volterra systems usfined for them (see in Hangos and Szederk´ e nyi (2012)) that ing a generalization of well-established methods available Volterra (LV) system. Dynamic similarity can also be destructure for a special subset of Lotka-Volterra systems usfined for them (see in Hangos and Szederk´ e nyi (2012)) that ing a generalization of well-established methods available fined for (see in Hangos and Szederk´ eenyi (2012)) that ing a generalization of well-established methods available shows thethem dynamic relationship between QP and chemical in the linear control theory. Of course, the results are not fined for them (see in Hangos and Szederk´ nyi (2012)) that ing a generalization of well-established methods available shows the the dynamic dynamic relationship relationship between between QP QP and and chemical chemical in control theory. Of the are shows in the the linear linear control theory. Ofa course, course, the results results are not not reaction 2009) with QP massand action law in supposed to be applicable forOf wide sub-class (within the shows thenetworks dynamic (Angeli, relationship between chemical the linear control theory. the results are not reaction networks (Angeli, 2009) with mass mass action action law supposed to be applicable for aa course, wide sub-class (within the reaction networks (Angeli, 2009) with law supposed to be applicable for wide sub-class (within kinetics innetworks a simple(Angeli, and transparent way.mass action law supposed Lotka-Volterra system class), therefore general(within results the are reaction 2009) with to be applicable for a wide sub-class the kinetics in a simple and transparent way. Lotka-Volterra system class), therefore general results are kinetics in aa simple and transparent way. Lotka-Volterra system general results are also given utilizing the class), notiontherefore of compartmental systems. kinetics in simple and transparent way. Lotka-Volterra system class), therefore general results are given utilizing the notion of compartmental systems. The control design of positive polynomial systems - to also alsothis given utilizing control the notion notion of compartmental compartmental systems. In generalized design case optimization based The control design of positive polynomial systems to also given utilizing the of systems. The control design of positive polynomial systems to generalized control which QP systems - has become quite popular The design belong of positive polynomial systems - to In In this thisare generalized control design design case case optimization optimization based based tools proposed. control whichcontrol QP systems systems belong has become quite quite popular In this generalized design case optimization based which QP belong --- has become popular tools are proposed. recently (see e.g. Tong et al. (2007)), thatquite is explained which QP systems belong has become popular tools are proposed. recently (see e.g. Tong et al. (2007)), that is explained tools are proposed. recently (see Tong et explained by the great practical andthat wideis expressive recently (see e.g. e.g. Tong importance et al. al. (2007)), (2007)), explained by the the great great practical importance andthat wideis expressive expressive by practical importance and wide byThis theresearch great is practical importance and wide expressive  partially supported by the Hungarian Research  This research is partially supported by the Hungarian Research  Thisthrough  Fund grant No. K-83440. research is supported Thisthrough research is partially partially supported by by the the Hungarian Hungarian Research Research Fund grant No. K-83440. Fund through through grant grant No. No. K-83440. K-83440. Fund

Copyright IFAC 2015 1010Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © 2015, IFAC (International Federation of Automatic Control) Copyright © IFAC 2015 1010 Copyright IFAC 2015 1010 Peer review© of International Federation of Automatic Copyright ©under IFAC responsibility 2015 1010Control. 10.1016/j.ifacol.2015.09.323

MICNON 2015 June 24-26, 2015. Saint Petersburg, Russia Attila Magyar et al. / IFAC-PapersOnLine 48-11 (2015) 1000–1005

2. LOTKA VOLTERRA MODELS AND THEIR UNDERLYING LINEAR DYNAMICS The so-called quasi-polynomial (QP) model is a set of nonlinear ODEs in the form   n m   Bik Aji zk , j = 1, . . . , m. (1) z˙j = zj Lj + i=1

k=1

where A ∈ Rm×n , B ∈ Rn×m are constant parameter matrices (coefficient matrix and exponent matrix, respectively), L ∈ Rm is a vector. The monomial-like terms in (1) of the form m  xi = zkBik , i = 1, . . . , n (2) k=1

are the so-called quasi-monomials of the system. It is known (e.g. Hern´ andez-Bermejo and Fair´en (1997)) that the set of quasi-polynomial systems can be split into classes of equivalence according to the matrix invariant M = B · A. (3) A unique, descriptive element of such QP equivalence class is the Lotka-Volterra model   m  (4) Mij xj  , i = 1, . . . , n x˙ i = xi λi + j=1

where M ∈ R , Λ = [λ1 , . . . , λn ]T ∈ Rn . Because of its descriptive nature, the values of the Lotka-Volterra parameter matrices can be computed from the QP parameter matrices of any quasi-polynomial model belonging to the same class of equivalence as M = B · A and Λ = B · L. The model (4) can also be written in a more compact matrixvector notation x˙ = diag(x1 , ..., xn ) M (x − x∗ ), (5) where diag(x1 , ..., xn ) is a diagonal matrix with xi in its ith diagonal element, x∗ is a unique positive equilibrium point of the system, which is the (nonzero) solution of the steady state version of (4): (6) 0 = Λ + M x∗ n×n

1001

It is important to note that the dynamic properties and the entire phase plot of the original and the transformed ODEs are the same including the steady state points and their stability properties, i.e. the time-reparametrization transformation is a similarity transformation. The matrices characterizing the time-reparametrized LotkaVolterra model are given as follows   ˆ = Im×m + 1(m+1)×1 · Ω (13) Aˆ = [ M Λ ] , B 01×m

where Im×m denotes the unit matrix of size m, and aj×k stands for the constant matrix of size j × k with elements a ∈ {0, 1}, and Ω = [Ω1 ... Ωn ] is the parameter vector of ˆ of the QP the transformation (7). The invariant matrix M model (13) is   M Λ ˆ =B ˆ Aˆ = ˆ M + 1(m+1)×1 · Ω · A. (14) 01×m 0 It is easy to see, that (14) depends linearly on the parameters Ωi of (7).

2.2 The translated X-factorable phase space transformation and the underlying linear dynamics Assume that the following set of ordinary differential equations (ODEs) dX = F (X) (15) dt is defined on the positive orthant P n . The singular solutions of Eq. (15) are defined by F (X) = 0. Consider the following nonlinear translated X-factorable transformation of Eq. (15) dX = diag(X1 , ..., Xn )F (X − C) (16) dt where the elements of C = [C1 , . . . , Cn ]T are positive real numbers, and X = [X1 , ..., Xn ]T .

Assume that F (X) is composed of polynomial-type functions with a finite number of singular solutions. It can be shown (Samardzija et al. (1989)) that the above transfor2.1 Time-reparametrization transformation mation can move the singular solutions into the positive There exists a nonlinear similarity transformation that orthant, and leaves the geometry of the state (or phase) preserve QP format, which is called time-reparameterization space unchanged within it (but not at or near the bound(Szederk´enyi et al. (2005)), time-rescaling or simply new- ary). Therefore, the dynamics of the solutions of Eqs. (15) time transformation (Hern´ andez-Bermejo et al. (1998)). It and (16) are structurally similar. introduces a nonlinear scaling of the time as follows. It is easy to see that a LV model has polynomial right-hand m  sides, so one can associate a structurally similar linear Ωk  zk dt (7) dt = ODE model k=1 (17) x˙ = M (x − x∗ ), Using (7) the original QP system (1) transforms to to the model (5), that is called the underlying linear n+1  m   dzj dynamic model of it. ˆ ˆji = z zkBik , j = 1, . . . , m. (8) A j dt i=1 k=1 2.3 Compartmental matrices where the number of quasi-monomials formally increase by ˆ ∈ R(n+1)×m as follows one, and thus Aˆ ∈ Rm×(n+1) and B Consider a linear ODE with constant coefficients where (9) the describing dynamic model is in the following form Aˆji = Aji , j = 1, . . . , m, i = 1, . . . , n Aˆj(n+1) = Lj , n = 1, . . . , m ˆji = Bji + Ωi , j = 1, . . . , n, B ˆ(n+1)i = Ωi , i = 1, . . . , m B

(10) i = 1, . . . , m (11) (12)

dx = Fx (18) dt Results from the theory of nonnegative and compartmental systems (see Chapter 2 of van den Hof (1996) and the

1011

MICNON 2015 1002 Attila Magyar et al. / IFAC-PapersOnLine 48-11 (2015) 1000–1005 June 24-26, 2015. Saint Petersburg, Russia

review paper Jacquez and Simon (1993)) will be used to analyze the stability properties of the above ODE. An n × n matrix F is called a compartmental matrix if it satisfies the following conditions

n  i=1

Fij ≥ 0, for i, j = 1, . . . , n, i = j Fij ≤ 0, for j = 1, . . . , n

(19)

Note that a compartmental matrix is a special type of Metzler matrix. The most important property of compartmental matrices is, that their eigenvalues are either zero, or they have strictly negative real parts. (In other words, compartmental matrices cannot have unstable or purely imaginary eigenvalues.) Moreover, Metzler matrices are known to be diagonally stable. Because F is also a Metzler matrix, it is also diagonally stable, i.e. there exists a positive diagonal matrix Q ∈ D+ such that ATk Q + QAk  0 (20) 2.4 Compartmental Lotka-Volterra models The key tool for forcing the sign pattern (19) for the Lotka-Volterra model having a positive equilibrium point (6) is time-reparametrization. As it was shown before, (7) is a similarity transformation that practically increases the chance of proving asymptotic stability. This is feasible when the Ω parameter vector of transformed LotkaVolterra model (14) can be selected in such a way, that ˆ the resulting Lotka-Volterra model coefficient matrix M admits the sign pattern (19). The conditions (19) can be derived for (14) as the linear constraints (21). n  Mii + Mji Ωj < 0, i = 1, . . . , n j=1 n 

Mij +

λj Ωj < 0,

j=1 n 

Mkj Ωk < 0,

k=1

λi +

n 

i = 1, . . . , n i = 1, . . . , n, j = 1, . . . , n, i = j

Because of the above strategy, the feedback is assumed to be a linear static state feedback u = u∗ − k (x − x∗ ) (23) 1×n ∗ where k ∈ R is a row vector, and x is the unique equilibrium point of the open loop system for a constant u∗ that satisfies (6). The closed loop system has the form   m m   Γi kj x∗j  (24) (Mij − Γi kj )xj + x˙ i = xi λ∗i + j=1

j=1

i = 1, . . . , n

3.1 Steady states and their invariance First we note that the equilibrium point of the closed loop system (24) coincides with that of the open loop system x∗ with the feedback (23). In vector notation, the steady state version of (24) is Using the same matrix-vector notation as before, the system (24) can expressed as follows. x˙ = diag(x1 , ..., xn ) M (x − x∗ ),

(25)

∗

where M = M − Γ k and x is the equilibrium point of the closed loop system. 4. CONTROLLER DESIGN The design is based on the fact, that nonlinear Xfactorable transformation described in subsection 2.2 enables us examine a dynamically similar linear ODE instead of the original one that can be characterized by the LTI matrix pair (M, Γ). 4.1 Single input case

(21) λj Ωj > 0,

i = 1, . . . , n

Mji Ωj > 0,

i = 1, . . . , n

j=1

n 

The motivation behind this simple structure is twofold: (i) the design will be based on the underlying linear dynamics that enables to have a linear static full state feedback for stabilizing, (ii) this structure corresponds to the model structure of lumped process models when the inlet flowrate is chosen as a manipulable input variable (see in (Hangos et al., 2004)).

j=1

Then a suitable time-reparametrization transformation Ω can be found by solving (21), if it exists. 3. LINEAR FEEDBACK STRUCTURE The (single) input extension of the system (4) is assumed to appear as follows   m  x˙ i = xi λi + Mij xj + Γi u , i = 1, . . . , n (22) j=1

This means xi u terms are appearing in the i-th state equation.

A necessary condition to be able to stabilize a dynamic system is its controllability (reachability). Because of the dynamic similarity of the original LV system and its underlying linear dynamics, it is enough to require that (M, Γ) is a controllable pair of a continuous time LTI system. This puts a structural condition of choosing the single input u for a given nonlinear QP system. Because of the special single input nature of the open loop LV system (22) and the chosen linear static feedback structure (23), one can in principle stabilize the system by using a simple linear pole-placement design (Astr¨ om and Wittenmark (1997)). In the sequel, the well-known linear quadratic (LQ) approach is used, i.e. the minimizing solution of   1 ∞ T J(x, u) = x Qx + uT Ru dt 2 0 (26) with respect to x˙ = M x + Γu

1012

MICNON 2015 June 24-26, 2015. Saint Petersburg, Russia Attila Magyar et al. / IFAC-PapersOnLine 48-11 (2015) 1000–1005

is computed via the continuous-time algebraic Riccati equation M T K + KM − KΓR−1 ΓT K + Q = 0 (27) −1 T and k = R Γ K is used in (23). The tuning parameters Q and R are supposed to be positive definite, symmetric matrices. This, however, will not be generalized for the multipleinput case, and may not guarantee diagonal stability, therefore a more general design is suggested below.

The first inequality stands for the non-negativity of the offdiagonal elements, i.e. (30), the second equality constraints describe the non-positivity of the column sums (31). Note, that there is still two questions to answer: How to enable negative kj and how to chose the parameter c of the objective function? The first problem can be solved using a trick generally used in LP problems: • For each unrestricted variable kj two variables are p n introduced: kj and kj • Everywhere in the LP the following substitution is used p n kj = kj − kj • The following nonnegativity constraints are used: p n kj ≥ 0 and kj ≥0 • In the solution of the modified LP p n · If kj = kj = 0, then kj = 0. p p n = 0, then kj = kj . · If kj > 0 and kj p n n > 0, then kj = −kj . · If kj = 0 and kj

4.2 General case The presented feedback structure (22) can be generalized to multiple input case and apply a more general controller design procedure. The multiple (p-dimensinoal) input extension is in the form: u = −k (x − x∗ ),  = 1, . . . , p (28) The closed loop system (4) is assumed to appear as follows:    p m   x˙ i = xi λi + Mij − Γi kj xj + j=1

+

p m   j=1 =1



Γi kj x∗j 

=1

i = 1, . . . , n

(29)

The primary aim of the design is to choose the feedback gain parameter matrix k ∈ Rp×m so that the closed loop matrix p  M =M− Γ k  =1

is compartmental (19), i.e. it has the following properties: n  i=1

M ij ≥ 0

i = j,

(30)

M ij ≤ 0,

j = 1, . . . , j

(31)

As (30-31) form a set of linear inequalities and (22) defines the k to appear linearly in M , a computationally feasible LP formulation is searched for. In order to apply a LP solver (Gass (1985)) for the above problem, first of all, it should be formulated as a problem in the form max cT y (32) such that A y ≤ b (33) and y ≥ 0 (34) where A, b, and c are constant matrices, vectors, respectively. Constraints The previous problem statement can be translated to the following inequality constraints, where the LP (decision) variable is the feedback gain, i.e. y is the vector consecutively containing the rows of the matrix k, i.e. y = kˆ = [k1 , . . . , kp ]. The linear constraints for the compartmental controller design are p  Γi kj ≤ Mij , ∀i = j, (35) =1

−

p  =1

kj

n  i=1

Γi ≤ −

n  i=1

Mij , ∀j = 1 . . . , n.

1003

This means, that the LP variable is the row vector kˆpn ∈ R2mp Objective function An obvious choice for the parameter c is the one that minimizes the control gain should be maximal so that the input energy applied should be minimal. On the other hand, this selection results in the most zeros in the feedback gains kj in the solution that ensures the minimal complexity of the resulting state feedback. In the constrained 2mp-dimensional parameter space it is   p n − kj + kj j,

that should be maximized, i.e. the chosen object function parameter is c = [−1, − 1, . . . , −1]T ∈ R2mp (37)

4.3 Extension with time-reparametrization A further extension can be introduced to the compartmental LV controller design problem. For the sake of simplicity, only the single input case is presented here. The joint problem of compartmental controller design (30-31) with the time-reparametrization based extension (21) introduces n extra parameters (Ω) to the problem and yields the LotkaVolterra coefficient matrix ˆ ∈ R(m+1)×(m+1) M

of the following form.

The main diagonal is in the form m ˆ = M − Γ k +  Ω (M − Γ k ), M ii ii i i j ji j i

i = 1, . . . , m

j=1

m

 ˆ M Ωj λ j (m+1),(m+1) = j=1

The off-diagonal elements are

(36)

m

 ˆ M Ω (Mj − Γ kk ), 1≤i=j≤m = Mij − Γi kj + =1

1013

MICNON 2015 1004 Attila Magyar et al. / IFAC-PapersOnLine 48-11 (2015) 1000–1005 June 24-26, 2015. Saint Petersburg, Russia

ˆ M (m+1),j =

m  =1

Ω (Mj − Γ kk ),

j = 1, . . . , m

m

 ˆ M Ω λ  , i,(m+1) =

i = 1, . . . , m

=1

The sign pattern (19) that ensures the compartmental property for the extended case is bilinear in the timereparametrization and the control gain parameters (Ω and k, respectively). The bilinear programming problem is a structured quadratic programming problem whose objective function is, in general, neither convex nor concave. The complexity of such problems can be NP-hard in the general case, i.e. they can be computationally demanding. 5. EXAMPLES 5.1 Nonlinear fermentation process example The example is a fermentation process system where the substrate and biomass can be fed to the reactor with constant concentrations SF and XF . The parameters of the system are listed in Table 1. The input of the model is the output flow-rate F . z˙1 = µmax z1 z2 + XF − F z1 (38) µmax z 1 z 2 + S F − F z2 z˙2 = − Y The model (38) is in QP form, it is apparent if it is written in the form (1):   z˙1 = z1 −F + µmax z2 + XF z1−1 (39)   µmax z˙2 = z2 −F − z1 + SF z2−1 Y In order to be apply the results of Section 4.1, the process model must be embedded into Lotka-Volterra form (for the details of the embedding procedure, see Hern´andezBermejo et al. (1998)).   µmax x˙ 1 = x1 −F − x 3 + SF x 4 Y x˙ 2 = x2 (F − µmax x1 − XF x2 )

2 1 1 1

The LQ matrices are selected in order that the closed loop system transients be faster:   10 0 0 0  0 10 0 0  Q= , R=1 (44) 0 0 10 0  0 0 0 10 The resulting linear static full state feedback gain is k = [−1.257 1.0526 − 2.5261 0.6879]. (45) The solution of the open loop and the closed loop dynamics can be seen in Figure 1, the corresponding inputs are given in Figure 2. It is easy to see, that the LQ controller has slightly increased the system dynamics due to the LQ parameters (44). 2.5

open loop closed loop

2 1.5 1 0.5 0

0.5

1

1.5 t

2

2.5

3

Fig. 1. Solutions of the open loop and the closed loop fermentation model (38). The dashed lines denote the open loop solution, the solid lines are the closed loop ones. It is apparent, that the controlled system is faster.

Table 1. Variables and parameters of the fermentation process model (38) biomass concentration substrate concentration output flow-rate substrate inlet feed flow-rate biomass inlet feed flow-rate yield coefficient kinetic parameter

Although (M, Γ) is not a fully controllable pair (the rank of the controllability matrix is 2), the uncontrollable eigenvalues of M are zeros according to the algebraic relationship between the state variables. Moreover, the fermentation model (38) is controllable. As the aim is to stabilize the (z1 , z2 ) dynamics which is according to (x3 , x1 ), the LQ based design can be used.

(40)

x˙ 3 = x3 (−F + µmax x1 + XF x2 )   µmax x˙ 4 = x4 F + x 3 − SF x 4 Y

z1 z2 F SF XF Y µmax ,

Linear state feedback The feedback structure used is the (23), i.e. it is in the form (42) F = F ∗ − k(x − x∗ ) The corresponding matrices of the Lotka-Volterra model are as follows.     0 0 −1 2 −1  −1 −1 0 0   1  M = , Γ= (43) 1 1 0 0  −1  0 0 1 −2 1

z

and

5.2 Example for the general case

[ gl ] [ gl ] [ hl ] g [ lh ] g [ lh ] 1 [h ]

The open loop QP system (39) has a locally stable positive equilibrium point with F ∗ = 1 at √ √ z ∗ = [ 2 − 1, 2 − 2]T ≈ [2.4142, 0.5858]T (41)

As an example for the generalized (LP based) multiple input controller design method suggested in Section 4.2, the following Lotka-Volterra system is regarded. x˙ 1 = x1 (5.385 − 3.859x1 + 0.879x2 + 1.049x3 + u1 ) x˙ 2 = x2 (−1.079 + 2.642x1 − 1.428x2 + 1.503x3 + u2 ) x˙ 3 = x3 (−3.11 + 2.067x1 + 2.163x2 − 1.927x3 + u1 ) (46) The system (46) has an unstable equilibrium point at

1014

u∗ = [0, 0]T ,

x∗ = [2.15, 1.87, 2.23]T .

(47)

MICNON 2015 June 24-26, 2015. Saint Petersburg, Russia Attila Magyar et al. / IFAC-PapersOnLine 48-11 (2015) 1000–1005

open loop input (F ∗ ) closed loop (control) input

F

1

0.5

0

0.5

1

1.5 t

2

2.5

3

Fig. 2. The input F of the fermentation model. The dashed line denotes the steady state input, the solid line is the control input. The applied linear feedback structure is   u1 u= = −k(x − x∗ ). u2

(48)

The following feasible solution has been obtained for the LP constraint set (35-36):   0.3967 0.120 0.09 k= (49) 0.2103 1.703 0.5167 Which means, that x∗ has been stabilized by means of a linear feedback via transforming the closed loop LotkaVolterra coefficient matrix to the compartmental matrix (50).   −4.256 0.759 0.959 M = 2.431 −3.131 0.987 (50) 1.670 2.043 −2.017 6. CONCLUSION A linear state feedback structure has been presented in this paper for a special subset of Lotka-Volterra systems. First, linear quadratic controller design method has been used for the stabilizing control of this nonlinear system class. Of course, the results is applicable only for a narrow sub-class (the set of Lotka-Volterra systems for which the corresponding QP models are controllable), therefore general results have also been given utilizing the notion of compartmental systems. In this generalized control design case optimization based tools are proposed. A further extension is also suggested that introduces extra degree of freedom via time-reparametrization, however, the problem to be solved turns to be a bilinear programming problem. ACKNOWLEDGEMENTS This research is partially supported by the Hungarian Research Fund through grant No. K-83440. REFERENCES Angeli, D. (2009). A tutorial on chemical network dynamics. European Journal of Control, 15, 398–406. Astr¨ om, K.J. and Wittenmark, B. (1997). Computercontrolled systems. Prentice Hall. Figueiredo, A., Gleria, I.M., and Rocha, T.M. (2000). Boundedness of solutions and Lyapunov functions in quasi-polynomial systems. Physics Letters A, 268, 335– 341.

1005

Gass, S.I. (1985). Linear Programming: Methods and Applications. McGraw-Hill. Hangos, K.M., Bokor, J., and Szederk´enyi, G. (2004). Analysis and Control of Nonlinear Process Systems. Springer-Verlag. Hangos, K.M. and Szederk´enyi, G. (2012). The underlying linear dynamics of some positive polynomial systems. Physics Letters A, 376(45), 3129 – 3134. Hern´andez-Bermejo, B. and Fair´en, V. (1995). Nonpolynomial vector fields under the Lotka-Volterra normal form. Physics Letters A, 206, 31–37. Hern´andez-Bermejo, B. and Fair´en, V. (1997). LotkaVolterra representation of general nonlinear systems. Math. Biosci., 140, 1–32. Hern´andez-Bermejo, B., Fair´en, V., and Brenig, L. (1998). Algebraic recasting of nonlinear ODEs into universal formats. J. Phys. A, Math. Gen., 31, 2415–2430. Isidori, A. (1995). Nonlinear Control Systems. SpringerVerlag. Jacquez, J. and Simon, C. (1993). Qualitative theory of compartmental systems. SIAM Review, 35(1), 43–79. Magyar, A. and Hangos, K. (2013). Control Lyapunov function based feedback design for quasi-polynomial systems. In 9th IFAC Symposium on Nonlinear Control Systems (NOLCOS 2013), Sep. 4-6, Toulouse. Magyar, A., Hangos, K., and Szederk´enyi, G. (2013). Stabilizing dynamic feedback design of quasi-polynomial systems using their underlying reduced linear dynamics. In 52nd IEEE Conference on Decision and Control (CDC 2013), Dec. 10-13, Florence. Magyar, A., Szederk´enyi, G., and Hangos, K.M. (2008). Globally stabilizing feedback control of process systems in generalized Lotka-Volterra form. Journal of Process Control, 18(1), 80–91. Samardzija, N., Greller, L.D., and Wassermann, E. (1989). Nonlinear chemical kinetic schemes derived from mechanical and electrical dynamical systems. Journal of Chemical Physics, 90 (4), 2296–2304. Szederk´enyi, G., Hangos, K., and Magyar, A. (2005). On the time-reparametrization of quasi-polynomial systems. Physics Letters A, 334, 288–294. Szederk´enyi, G., Lipt´ak, G., Rudan, J., and Hangos, K. (2013). Optimization-based design of kinetic feedbacks for nonnegative polynomial systems. In IEEE 9th International Conference of Computational Cybernetics, July 8-10, Tihany, Hungary, 67–72. doi: 10.1109/ICCCyb.2013.6617563. ISBN: 978-1-47990063-3. Tong, C.F., Zhang, H., and Sun, Y.X. (2007). Controller design for polynomial nonlinear systems with affine uncertain parameters. Acta Automatica Sinica, 33(12), 1321 – 1325. doi:10.1360/aas-007-1321. van den Hof, J.M. (1996). System Theory and System Identification of Compartmental Systems. Ph.D. thesis, University of Groningen. VanAntwerp, J. and Braatz, R. (2000). A tutorial on linear and bilinear matrix inequalities. Journal of Process Control, 10, 363–385.

1015