Home

PnPMPC Toolbox v. 1.0 - User manual

1. 24 e x0 is the initial local state x 0 a vector of dimension n x 1 e din is the vector of exogenous inputs which act on the subsystem and they are not distur bances Indices of these exogenous inputs are in ob jcenmpc Exo e xrefin means the reference state of nominal system over the control horizon it can be expressed as a vector n x 1 reference constant over the control horizon or a matrix n x N the columns 1 refer to the values at the time instant 1 within the control horizon columns 2 refer to the values at time 2 etc If this parameter is empty or not passed it is set to zero by default e urefinin means the reference input of nominal system over the control horizon it can be expressed as a vector m x 1 reference constant over the control horizon or a matrix m x N the columns 1 refer to the values at the time instant 1 within the control horizon columns 2 refer to the values at time 2 etc If this parameter is empty or not passed it is setted to zero by default e options sdpsettings object for YALMIP options As output we have e control input value e diagnostics information about the optimization problems 4 2 Examples Examples of cenmpc class can be found in the examples cenmpc directory Execute examplelcenmpc for a deterministic MPC controller and example2cenmpc for tube based MPC 25 Chapter 5 The pnpmpc class 5 1 Introduction In this chapter we introduce the pnpmpc class i
2. Moreover ujj and uf are forces applied to the first and second mass respectively We also have Au Ao A21 Ago 0 1 A Ay Ag k h M M M M 0o 0 Ar z x a 3 2 B diag B11 B22 B Bo 4 3 3 M First of all we present the methods that will be used 11 The lss method lss is the constructor of the homonymous class with this declarations objlss 1ss A B ni mi C D pi name Ts where A and B are the matrices in 2 7 ni and mi are vectors storing number of local states and local control inputs for each subsystem Cand D are the matrices in 2 8 pi is a vector storing number of outputs for each subsystem name is an optional string defining the system name Ts is the sampling time that can be Ts for continuous time system Ts gt 0 for discretized system with sampling time Ts Ts 1 for system already in the discrete time form It is possible to introduce directly the matrices of the collective system system 2 7 and 2 8 but this operation is usually time consuming and error prone Hence in the example we will start from an empty object and add subsystems with addSystem method The addSystem method The addSystem method adds to an 1ss object a single subsystem The method declaration is objlss objlss addSystem Ai Bi Ci Di name Ts where the subsystem matrices are Ai Bi Ci Di as defined in 2 5 and 2 6 The arguments name and Ts
3. npoints 500 X zonotope zeros 2 1 2 0 0 3 xf f Ox wi CL zeros dimension X npoints 1 0 1 x 1 0 5 x 2 exp 0 1 x 1 72 0 1 0 9 x 1 0 1 x 2 0 1 cos x 2 0 05 x 2 72 for i 1 size xf 2 xx randpoint X xf i end definition of g and h such that f g h where g and h are convex functions 0 1 0 9 x 1 0 1 cos x 2 0 05 x 2 72 J x w L x w L g h options for options split options split options split compute zonotope approximation of f X f xx 0 5 x 2 1 0 1 x 1 exp 0 1 x 1 72 outerApproximation function ngenerators 2 alpha 0 5 max_zono 5 68 set X and sampling of f X 0 1 x 2 1 1 Zu approximation using 1 zonotope ZuVec approximation using options split max_zono zonotopes Xs zonotope X splitted J jacobian and H hessian are outputs variables for future computations Zu ZuVec Xs J H outerApproximation X f g h options 4 plots figh figure 1 subplot 1 2 2 hold on Zu plotZ y ZuVec plotZ r plot Cat 1 1 xf 252 5 Vba box on subplot 1 2 1 Xs plot g box on 1 w N I N Figure 12 6 Results of example example2outerA pproximation m 12 3 3 example3outerA pproximation m definition of function f set X and sampling of f X differently from example2outerApproximation we want to split only along x
4. 9m are not linearly independent Given G e R with Rank G n then the bisection is complete i e Bisect Z provides two sub zonotopes that do not overlap The operators Bisecty and Splitx have been implemented in PnPMPC toolbox in the following two instructions 62 Z bisect Z split gener alpha The following code shows an example where zonotope Z is splitted P G alpha 2 zonotope p G Z bisect Z split k alpha figure 1 Z plotZ figure 2 Znb plotZ struct shade 0 6 figure 3 Zns plotZ struct shade 0 6 Figure 12 2 An example of split Bisect Z 63 12 2 2 Bounding box operator The smallest interval vector containing Z and having its same center is given by the operator rs row sum Given a matrix G R rs G is a n x n diagonal matrix m rs G ii 5 Gi l Liu yh j l The operator rs has been implemented in PnPMPC toolbox in the following instruction BB Z bounding_box where BB is a zonotope object computed using operator rs The following code shows an example where zonotope Z is bounded by the smallest interval vector do y Is 1242 45 31456 2451 5 zonotope p G Z bounding_box figure 1 hold on BB plotZ Z plotZ b Figure 12 3 An example of bounding box of a set Z 64 12 2 3 Reduction operator Another important operator is the reduction operator whose purpose is to outer bound a given zonotope wit
5. Mi the control scheme is completely decentralized since uj depends on local quantities only The key condition enabling PnP design is the following a gt Si Ae FBK Ay 6 BKig lloo lt 1 JEN k 0 Therefore in order to design an unconstrained PnP controller for subsystem we need to solve the following problem Problem P Given 6 compute matrices Ky and Kij j NM such that a lt 1 If Problem P cannot be solved we declare that it is impossible to design an unconstrained PnP controller for subsystem The interested reader is refer to 19 and Chapter 4 in 4 7 1 1 The pnpctrl method Assume subss stores the model of subsystem i given by 2 3 and 2 4 Controller Cj is created through 43 ctrl pnpctrl subss Xj flag The input argument flag is a structure with the following possible fields e distributedParents indices of subsystems such that 6 1 in 7 1 It must be a subset of the indices of parent subsystems e norm specify norm for the computation of matrices K similar as in Chapter 7 in 4 boundKi j specify bound for the computation of matrices K similar as in Chapter 7 in 4 Qi and Ri Qi and Ri are inputs for the nonlinear optimization problem as proposed in Chapter 4 and 5 in 4 e useConstraints if false the optimization problem does not take advantages of the knowledge of scaling factors for states of parent subsystems t rue otherwise Scaling factors ar
6. and w t e W C R ifVa t X there exist u t U such that x t 1 X Vw t e W 10 1 e mRPI In this section we explain how to use the epsilon_mRPI class in order to compute an e outer approximation of the mRPI set for a linear constrained systems Our implementation is based on the algorithm proposed in 27 We propose two examples 10 1 1 examplelepsilon_mRPI m The code below shows how to compute an e mRPI set for the LTI system at Ax Bu w where x R u Kz and w lies in the polytope W Non constraint on the state is assumed and hence X R matrices of the LTI system A 11 01 53 wW Il 145 1 1 1 17 1 03 Il bounded disturbances as Polyhedron object Polyhedron eye 2 eye 2 1111 approximation epsilon 5 10 5 create an object F that is a epsilon mRPI for the closed loop LTI system with dynamic A BK and disturbances bounded in W W is a Polyhedron object epsilon_mRPI A B K W epsilon Otherwise one can create an object F that is a epsilon mRPI for the closed loop LTI system with dynamic A BK and disturbances bounded in W W is a zonotope object W zonotope 0 5 0 1 eye 2 F epsilon_mRPI A B K W epsilon plot the approximation of the mRPI plot Other useful methods fot the epsilon _mRPI class are isinside Test if a point is inside the e mRPI set F E dist Js F isinside x doubleHK If W is a Polyhedron object th
7. all controllers will be designed 45 Chapter 8 The pnpEstimator class 8 1 Introduction In this chapter we introduce the pnpEstimator class in order to design a PnP state estimator for a single subsystem PnP state estimators are described in 6 and Chapter 8 in 4 In those papers the authors propose a distributed for 2 5 and 2 6 We assume that local outputs depend only on local states and disturbances We define for i M the Local State Estimator LSE as Sa Zh Andy Baug Lalu Cit DI Agee DI dala Gi GEN JEN Se Biju BcenUcen Mid JENG 8 1 where j R is the state estimate L i R are gain matrices and Oi 0 1 This implies that Zi depends only on local variables j 4j uj and yj and parents variables fiz ujj Yy J E Ni centralized inputs ucen and exogenous inputs d Binary parameters big ZEN can be chosen equal to one for exploiting the knowledge of parents outputs or equal to zero for reducing the number of transmitted output samples The key condition enabling PnP design is the following a DI DO IIA Lat Ay big LigCy loo lt 1 JEN k 0 Therefore in order to design a PnP LSE for subsystem i we need to solve the following problem Problem P Given dij compute matrices Li and Lij j NM such that a lt 1 If Problem P cannot be solved we declare that it is impossible to design a PnP LSE for subsystem 2 The interested reader is refer
8. 3 1 Memory optimization Taking into account that usually LSS have a sparse structure and therefore matrices in 2 7 and 2 8 have many zero elements we decided to use the sparse matrices of MatLab by default This usually allows one to save a considerable amount of memory To visualize the matrices in the full format the MatLab command full can be used To save further memory we also used the following trick If in 2 8 one has y x then e the matrix C is the identity matrix of order n e Dis a matrix of zeros R e there are no exogenous inputs or centralized control inputs acting on the output Therefore we decided to not save these matrices However if the user requires them with a get method they will be generated upon demand In the following this particular form of the output equation will be called the standard setting Model discretization Continuous time models can be discretized using the method clss2d1ss that means continuous LSS to discrete LSS with this declaration dobjlss objlss clss2dlss Ts method where Ts is the sampling time and method is the discretization technique All the discretiza tion techniques of the function c2d of the Control System Toolbox which converts continuous time systems to discrete time models are implemented like zoh zero order hold or Tustin see help c2d for more information The c2d discretization methods are not always a good choice for LSS because for exact d
9. Iss object are created with the method createPnPEstimators described next 8 2 1 The createPnPEstimators method This method acts on an Iss object and it is called as follows estimators L createPnPEstimators objlss flag where flag has the same meaning as in the pnpEstimator method This method extracts a subsystem calls the constructor of class pnpEstimator and design the state estimator The output estimators is an array of dimension 1xnumSys and in the i th position the state estimator for subsystem i is stored L is the overall matrix for the LSS collecting matrices Lij 8 3 Examples In this section we propose some examples Section titles corresponds to m files that can be found in the examples pnpEstimator directory 8 3 1 examplelpnpestimator m Executing file examplelpnpestimator the example proposed in 6 and Chapter 8 in 4 is generated designing the LSS the state estimators and simulating the convergence of the state estimation Note that since the LSS as several states for each subsystem the generation of state estimators could take minutes 49 Chapter 9 The lse class 9 1 Introduction In this chapter we introduce the lse class in order to design a large scale state estimator for an LSS PnP state estimators are described in 5 and Chapter 3 in 4 In those papers the authors propose a distributed for 2 5 and 2 6 We assume that lo cal outputs depend only on local s
10. The reduction operator redngen Z is applied to the green zonotope yielding a more conservative approximation red zonotope 12 2 5 Zonotope Polyhedron e Polyhedron P is a zonotope bool iszonotope P e Polyhedron P to Zonotope Z Z Polyhedron2zonotope P e Zonotope Z to Polyhedron P P zonotope2Polyhedron Z 12 3 Outer approximation of a nonlinear function The function outerApproximation allows one to compute an outer approximation of a non linear function evaluated on a zonotope Z We compute the outer approximation of the nonlinear function using the methods proposed in 31 30 In 31 an algorithm to compute zonotope outer approximations for nonlinear systems was proposed The authors suggest to create an image of a zonotope through a nonlinear function using DC programming which is based on DC functions 1The function outerApproximation needs Symbolic Math Toolbox 66 A DC function f x R R is a function that can be expressed as the difference of two convex functions i e f x g a h x where g x and h x are convex functions In order to compute a tighter outer approximation in 30 and 32 the authors proposed an algorithm in order to split the zonotope set Z where the function f x is more nonlinear and then the outer approximation is obtained as the union of the outer approximations obtained using the DC programming on each zonotope of the splitting In the following we propos
11. application that can benefit of decentralized and distributed control schemes is the regulation of a Power Network System PNS Here we describe the PNS proposed as a benchmark exercise 3 within the HYCON2 project 23 We consider a PNS as composed by several power generation areas coupled through tie lines 24 The aim is to design the Automatic Generation Control AGC layer for frequency control with the goal of e keeping the frequency approximately at the nominal value e controlling the tie line powers in order to reduce power exchanges between areas In the asymptotic regime each area should compensate for local load steps and produce the required power We consider thermal power stations with single stage turbines The dynamics of an area equipped with primary control and linearized around equilibrium value for all variables can be described by the following continuous time LTI model 24 of i Aura Bux LAP 5 Aij g 6 1 GEN where x A0 Awi APmi AP is the state uj APref is the control input of each area AP is the local power load and N is the sets of neighbouring areas i e areas directly connected 34 to oh through tie lines The matrices of system 6 1 are defined as 0 1 0 0 0 Dien Pis _ Di 1 0 0 Aii Pij jen Hi B III 0 0 7 Ti 0 0 0 7 RiTg Tg Toi 6 2 0 000 0 Pij Aij 2H 0 0 0 L mn IF 0 000 0 0 000 0 For the meaning of constants as
12. control invariant set plot Other useful methods for the parameterizedRCI class are 56 isinside Test if a point is inside the parameterized robust control invariant set R 1 1 R isinside x ulnv Given x X compute the control action u x such that x t 1 X Vw W zeros 2 1 R uInv x double The function returns the H representation of the parameterized robust control invariant set Moreover a Polyhedron object of the set can be obtained through the attribute RP R double R RP 10 4 Zonotopic robust control invariant sets In this section we explain how to use the zonotopeRCI class in order to compute a zonotopic robust control invariant set for a linear constrained systems Our implementation is based on the algorithm proposed in 28 We illustrate the use of the class zonotopeRCI in the following example 10 4 1 examplelzonotopeRCI m matrices of the LTI system 11 01 1 0 11 matrices for constraints on state and input variables X x cx x lt dx u cu u lt du O 1 0 0 a s 285 3 3 8 42 2 Ja y 1 Ss 3 3 parameter of the algorithm greater than the number of states 5 matrices of the zonotope for the description of set W w f Ed d inf lt E 10 011 57 f O i 0 l create the zonotopic robust control invariant set qa 1 qb 1 gp i Z zonotopeRCI A B cx dx cu du f E k qa qb qp 4 plot the zo
13. ejj J N are in the set Ejl ejl Hej1e71 lt Kej Ej2 ej2 H 2ej2 Kej2 where ejl and ej2 are state estimation errors of two parent subsystems of the local subsystem subss These sets are needed if flag boundedEstimation true e Cj is a cell matrix C1 Cj2 where each element is a the output linear transformation of parent subsystems Matrix Cj is needed if di 1 e Oj is a cell matrix Hoj1 Koj1 Hoj2 Koj2 such that parent disturbances dij j N acting on outputs of parent subsystems are in the set Oj1 dj1 Hoj1djl lt Koji 0j2 dj2 Hoj2dj2 lt Koj2 where djl and dj2 are disturbances of two parent subsystems of the local subsystem subss These sets are needed if flag boundedEstimation true The input argument flag is a structure with the following possible fields e distributedParents indices of subsystems such that d 1 in 8 1 It must be a subset of the indices of parent subsystems e norm specify norm for the computation of matrices L similar as in 6 and Chapter 8 in 4 e Qi and Ri Qi and Ri are inputs for the nonlinear optimization problem as proposed in Chapter 8 in 4 47 e boundedEstimation if false the optimization problem design a state estimator guar anteeing bounded error state estimation t rue otherwise e ExoBounded indices of exogenous inputs acting on subsystem subss that are distur bances therefore if boundedEstimation true pnp
14. for efficient robust model predic tive control of constrained linear discrete time systems subject to bounded disturbances in Proceedings of the 16th IFAC World Congress Prague Czech Republic July 4 8 2005 pp 241 246 S Riverso M Farina and G Ferrari Trecate Plug and Play decentralized Model Predictive Control in Proceedings of the 51st IEEE Conference on Decision and Control Maui Hawaii USA December 10 13 Dec 2012 pp 4193 4198 S Riverso and G Ferrari Trecate Plug and Play distributed model predictive control with coupling attenuation Optimal Control Applications and Methods p Submitted 2013 S V Rakovi and M Baric Parameterized Robust Control Invariant Sets for Linear Sys tems Theoretical Advances and Computational Remarks IEEE Transactions on Automatic Control vol 55 no 7 pp 1599 1614 2010 S Boyd L El Ghaoui E Feron and V Balakrishnan Linear matrix inequalities in system and control theory Philadelphia Pennsylvania USA SIAM Studies in Applied Mathematics vol 15 1994 IBM IBM ILOG CPLEX Optimization Studio 12 4 2011 Hycon2 Highly complex and networked control systems HYCON2 Network of excellence 2010 Online Available http www hycon2 eu H Saadat Power System Analysis 2nd ed New York NY USA McGraw Hill Series in Electrical and Computer Engineering 2002 S Riverso M Farina and G Ferrari Trecate Design of plug and play
15. lt Kaj 2 23 and for all ucen jj j E 1 Nnu it holds Veen ucen j E Rs Heen jUcen j lt Keen 2 24 Coupling constraints have a totally analogous representation Chapter 3 Modeling of LSS the Iss and subss classes 3 1 The Iss class There are many MatLab toolboxes that offer facilities for modeling MIMO and LTI systems However most of them have not been designed to handle the interconnection of several subsystems Therefore it is not immediate to manage an LSS that means to add remove subsystems to extract a subsystem to set coupling terms etc In the PnPMPC toolbox we have implemented methods for these tasks so as to ease the modeling of systems in the form 2 7 and 2 8 Moreover the system object stores the constraints and some other useful pieces of information The main properties of an 1ss object are numSys the number s of subsystems composing the LSS Ts the sampling time all matrices appearing in 2 7 and 2 8 ni mi pi vectors of s elements which collect the numbers of states inputs and outputs of every subsystem as in 2 9 2 10 and 2 11 numExo numICen are respectively nq in 2 12 and nu 2 13 all constraint matrices H and K appearing in 2 21 Hdeltau Kdeltau that will be explained in Section 3 2 4 couplingmatrix name that will be explained in Section 3 1 namex nameu namey named nameucen in order to assign names to variables see Section
16. model predictive con trol an approach based on linear programming in Proceedings of the 52nd IEEE Conference on Decision and Control Florence Italy December 10 13 2013 pp 6530 6535 S Dashkovskiy B S Riiffer and F R Wirth An ISS small gain theorem for general networks Mathematics of Control Signals and Systems vol 19 no 2 pp 93 122 2007 S V Rakovi E C Kerrigan K I Kouramas and D Q Mayne Invariant approximations of the minimal robust positively invariant set IEEE Transactions on Automatic Control vol 50 no 3 pp 406 410 2005 73 28 29 30 31 82 S V Rakovi Robust Control of Constrained Discrete Time Systems Characterization and Implementation Ph D dissertation Imperial College London University of London 2005 D Barcelli N Bauer and P Trnka WIDE Toolbox 2012 Online Available http ist wide dii unisi it index php p toolboxsp D M Raimondo S Riverso S Summers C N Jones J Lygeros and M Morari A set theoretic method for verifying feasibility of a fast explicit nonlinear Model Predictive Controller in Distributed Decision Making and Control R Johansson and A Rantzer Eds Springer Lecture Notes in Control and Information Sciences vol 417 2012 ch 13 pp 289 311 T Alamo J M Bravo M J Redondo and E F Camacho A set membership state es timation algorithm based on DC programming Au
17. of the PnPMPC toolbox for modeling purposes will be illustrate in Chapter 3 2 1 Large scale linear systems In this toolbox we consider Multi Input Multi Output MIMO Linear Time Invariant LTT systems either in discrete or continuous time An LSS is composed by s physically interconnected subsystems In the following dependence of variables from time will be omitted except where indicated The dynamics of subsystem i M 1 s is ct Aux Baug DO Ag Big 2 1 SEN where zp R is the state at time k zti is the state at time k 1 for continuous time systems Ti stands for Tta and ur R is the local input at time k Moreover we have A R Ai R i By R and By R Next we define the index set M appearing in 2 1 Definition 3 Subsystem i is dynamically coupled to subsystem j if Aij 0 Definition 4 Subsystem i is input coupled to subsystem j if Bij A 0 It is possible to use graph theory to represent the coupling between the subsystems The coupling graph of a system is directed graph where nodes are the subsystems and the set of edges is given by i j i j Aij 7 0 or Bij 0 An example is given in Figure 2 1 Coupling allows us to define the predecessors of a subsystem Definition 5 The set of parents to subsystem i is N j i j x 1 SubSystem 1 SubSystem 3 x 3 x 2 SubSystem 2 SubSystem 4
18. of the continuos linear system pnsCn where n is the number of the scenario e parameters of the control experiment Tsim and aPL where aPL corresponds to AP e the results of the control experiment x aPref eta and Phi where aPref corresponds to APref For each Scenario we included also a Simulink model In particular one can load the file mat of a control experiment and simulate the power network system given the power load steps and the power reference computed through centralized MPC or PnPMPC controllers 41 6 4 1 Example of simulation In the following we illustrate how to use the files mat and the Simulink models for designing centralized and pnpmpc controllers e Step 1 We can simulate different scenarios using the Simulink models present in the folder of each scenario For Scenario 2 we then open the file simulatorPNS_AGC_2 mdl This step is performed with the MatLab command open simulatorPNS_AGC_2 e Step 2 Centralized case In the subfolders of each scenario there are MatLab files for centralized simulations as scenario number scenario D Dss Diag Full Zero Data Assume we want to simulate Scenario 2 using the discretization Dss and centralized MPC with zero terminal constraint M PCzero We need to load MatLab file scenario2DssZeroData This operation can be performed with the MatLab command load scenario2DssZeroData PnPMPC case In the subfolders of each scenario there are MatLab files for PnPMPC sim
19. similar way In addition many properties are the same as those of the lss class Next we describe only the new ones referring the reader to the html documentation in the doc directory for a comprehensive description of all methods and properties If we extract subsystem i from an Iss object we lose coupling with other subsystems Therefore we created the following properties to save e couplingA couplingB couplingC couplingD are scalar vectors containing the in dices of the subsystems in M that are coupled with the extracted subsystem i e subsystem i couplingA will contain the indices of the subsystems coupled through the A matrix couplingB couplingC and couplingD have a similar meaning e Aij Bij Cij Dij Aij contains blocks A where i is the number of the extracted subsystem and 7 couplingA Same remarks apply to the other properties 21 3 3 1 examplelsubss m In this example we show how to add a subsystem to an lss object using a subss object example2lss A 0 1 k M h M B 0 100 M sub subss A B sub sub addStateConstraint 3 5 0 1 1 2 modelCart modelCart addSystem sub 3 3 2 example2subss After running example3lss we can extract subsystem 237 as follows example3lss sub modelCart getSystem 237 With the previous command subsystem sub is an object of the class subss In this case only local constraints are saved For example constraints among subsystem 237 and 238 wil
20. the scenarios introduced in Section 6 1 Different control schemes will be compared with the centralized MPC scheme described next For a given Scenario at time t we solve the centralized optimization problem P x t 6 3a t N 1 elit DI x la 8 ulle N x s 6 36 x k 1 Ax k Bu k LAPrL t ke0 N 1 6 3c x k X ke0 N 1 6 34 u k U kKe0 N 1 6 3e x N e Xf 6 3f and then apply AP e u 0 We note that the cost function depend upon x and u that are defined as xh 0 0 APr APr and us APr The constraints X and U in 6 3d and 6 3e are obtained from constraints listed in Table 6 2 In the cost function 6 3b we set N 15 Q diag Q1 Qar and R diag R1 Rm where 500 0 0 0 0 001 0 0 Qi and R 10 0 0 001 0 0 0 0 10 Weights Q and R have been chosen in order to penalize the angular displacement A0 and to penalize slow reactions to power load steps Since the power transfer between areas i and j is 38 given by APrie k Pi A0 k A0 k 6 4 the first requirement also penalizes huge power transfers In order to guarantee the stability of the closed loop system we design the matrix S and the terminal constraint set Xy in three different ways e S is full MPC full we compute the symmetric positive definite matrix S and the static state feedback auxiliary control law Kauxx by maximizing the volume of the ellipsoid de scribed by S insi
21. vol 58 no 10 pp 2608 2614 2013 Plug and Play Model Predictive Control based on robust control invariant sets Au tomatica p Accepted 2014 Plug and Play Model Predictive Control based on robust control invariant sets Dipartimento di Ingegneria Industriale e dell Informazione Universita degli Studi di Pavia Pavia Italy Tech Rep 2012 Online Available arXiv 1210 6927 S Riverso Distributed and plug and play control for constrained systems Ph D disserta tion Universita degli studi di Pavia 2014 S Riverso D Rubini and G Ferrari Trecate Distributed bounded error state estimation for partitioned systems based on practical robust positive invariance in Proceedings of the 12th European Control Conference Zurich Switzerland July 17 19 2013 pp 2633 2638 S Riverso M Farina R Scattolini and G Ferrari Trecate Plug and play distributed state estimation for linear systems in Proceedings of the 52nd IEEE Conference on Decision and Control Florence Italy December 10 13 2013 pp 4889 4894 M Herceg M Kvasnica C N Jones and M Morari Multi Parametric Toolbox 3 0 in Proceedings of the 12th European Control Conference Zurich Switzerland July 17 19 Jul 2013 pp 502 510 Online Available http control ee ethz ch mpt J L fberg YALMIP A toolbox for modeling and optimization in MATLAB in Proceedings of IEEE Symposium on C
22. well as some typical parameter values we refer the reader to Table 6 1 Deviation of the angular displacement of the rotor with respect to the stationary reference axis on the stator Speed deviation of rotating mass from nominal value Deviation of the mechanical power from nominal value p u Deviation of the steam valve position from nominal value p u Deviation of the reference set power from nominal value p u Deviation of the nonfrequency sensitive load change from nominal value p u Inertia constant defined as H Kinetic energy at rated epee typically values in range 1 10 sec machine rating Speed regulation Defined as Percent change in load change in frequency Prime mover time constant typically values in range 0 2 2 sec Governor time constant typically values in range 0 1 0 6 sec Slope of the power angle curve at the initial operating angle between area i and area j Table 6 1 Variables of a generation area with typical value ranges 24 p u stands for per unit We note that model 6 1 is input decoupled since both AP and APz act only on subsystem DE Moreover subsystems Di are parameter dependent since the local dynamics depends on the Vjen Pij 2H i In the following we introduce three scenarios corresponding to different interconnection topolo quantities gies of generation areas The model parameters and constraints on A and on APref for systems in all Scenar
23. with the following constraints cy EX CR 2 162 ujj U CR 2 16b We furthermore assume that X U and Y are polyhedra It is also possible to define coupling constraints between two or more subsystems In the case of two subsystems and j we have z cp Xiz CRO 2 17a up Uy Uiz C RT 2 17b Yip n Yiz E RTP 2 17c la Yea Yay Constraints for the global lss system can be defined as follow x XCR 2 18a ue UCR 2 18b yeYCRP 2 18c where X I X 0 Q Xy 2 19a iEM i jEM i j U U N Q Uy 2 19b ieM i jEM i j Y I ZA Q Yo 2 19c ieM i jEM i j Moreover we can define constraints for the exogenous inputs and for the centralized control inputs deDCR 2 20a Ucen Ucen C Ru 2 20b Recalling Definition 1 we have matrices Hy Hu Hu Ha Heen and vectors Kz Ky Ku Ka Keen of suitable dimensions such that X x ER H x lt Ky 2 21a U u R Hu lt Ku 2 21b Y y R Hy lt Ky 2 21c d e R 4 Had lt Ka 2 21d Ucen Ucen ER HoenUcen lt Keen 2 21e Moreover for each subsystem i M there are matrices Hx Hu Hy and vectors Kz Ku Ky such that X zq ER Ay 2 lt Ka 2 22a U uy E R Hug lt Ku 2 22b Yi yp E R Ay yg lt Kyi 2 22c In a similar way for dj 7 1 na we have j dij ER Ha jdij
24. x 4 n n Figure 2 1 Example of a coupling graph of a system composed by four subsystems A black arrow from system to system j means that A 0 a red arrow means By 0 For example from Figure 2 1 one has M2 1 4 since A21 B21 and Agq are not null matrices We also define the set of children of subsystem i e the set of subsystems influenced by it Definition 6 The set of children to subsystem i is Si j j i E For example from Figure 2 1 one has S2 4 since A42 4 0 The output of system is given by Yi Curly Diuji 5 Citi Diju 2 2 SEN where yj RP Cu RP Dy RE Ci RP and Di e REX Definition 7 Two subsystems i and j are output coupled if Ci or Dij are different from zero It is also possible to include exogenous signals that act on system i by replacing 2 1 and 2 2 with JEN IENA vi Cic Daug DO Cie Dig DI Nidy 2 4 GEN JENa where djy R is an exogenous input Na C N is the set of exogenous inputs that act on subsys tem i and Mij R Nij RP We further enhance our model by introducing centralized control inputs cen 80 that the dy namics of subsystem 7 becomes zti Aut Buug DO Ac Bgu XO Myd XO Boenijtcenj 2 5 SEN IENAa i SEN ui yt Cisti Daug DO Cig Disp XO Nady DI DoensisUcenta 2 6 SEN JEN di JEN uri where N i C N is the set of centralized control inputs Ucen j J R t
25. 1 therefore x 2 is substituted by w 1 npoints 500 X1 zonotope 0 2 69 zonotope 0 3 X X1XW xf zeros dimension X npoints f x w 1 0 1 x 1 0 5 w 1 exp 0 1 x 1 72 0 1 0 9 x 1 0 1 w 1 0 1 cos w 1 0 05 w 1 72 for i 1 size xf 2 xx randpoint X xf i f xx 1 xx 2 end 4 definition of g and h such that f g h where g and h are convex functions g x w 0 5Xw 1 0 1 0 9 x 1 0 1 cos w 1 0 05 w 1 72 1 h x w L 1 0 1 x 1 exp 0 1 x 1 72 O 1 w 1 1 options for outerApproximation function options split ngenerators 2 options split alpha 0 5 options split max_zono 5 compute zonotope approximation of f X Zu approximation using 1 zonotope 4 ZuVec approximation using options split max_zono zonotopes Xs zonotope X splitted J jacobian and H hessian are outputs variables for future computations Zu ZuVec Xs J H outerApproximation X1 f W g h options 4 plots h figure 1 subplot 1 2 2 hold on Zu plotZ y ZuVec plotZ r plot xt ist 5 xt 252 5 Ba box on subplot 1 2 1 Xs plotZ g box on 70 N Figure 12 7 Results of example example3outerApproximation m 71 Bibliography 1 CI 10 11 12 S Riverso M Farina and G Ferrari Trecate Plug and Play Decentralized Model Predictive Control for Linear Systems IEEE Transactions on Automatic Control
26. Estimator design a robust state estimator as in 6 and Chapter 8 in 4 estimator is the PnP state estimator designed with options given in flag 8 1 3 The computeEstimation method The state estimation jj is computed as in 8 1 executing the function xest yest yfilt estimator computeEstimation xt yt ut xtj ytj utj utcen d e xt is Zj at time t e yt is yj at time t e ut is uj at time t e xtj is a vector of x1 7 N at time t e ytj is a vector of yj j E N at time t needed if bi e utj is a vector of uj j N at time t needed if subsystems i and j are input decoupled e ucen iS Ucen at time t e dis dat time t i e exogenous inputs that does not act on the subsystem as not measurable disturbances As output we have e xest is in at time t 1 e yest is yj estimated at time t e yfilt is yj estimated at time 1 8 1 4 The checkErrorEstimation method If flag boundedEstimation true this method checks if a given initial error ejj 0 E and hence if the state estimator Zi can guarantee convergence of the state estimation and bounded error at all time instants test estimator checkErrorEstimation e 8 1 5 The pnpEstimator2ss method Given a PnP state estimator returns a state space object sys estimator pnpEstimator2ss sys is an ss object 48 8 2 Design of PnP state estimators for a large scale system PnP state estimators Zi for each subsystem in an
27. PnPMPC Toolbox v 1 0 User manual Stefano Riverso Alberto Battocchio and Giancarlo Ferrari Trecate Dipartimento di Ingegneria Industriale e dell Informazione Universit degli Studi di Pavia via Ferrata 1 27100 Pavia Italy April 2014 The research leading to these results has received funding from the European Union Seventh Frame work Programme FP7 2007 2013 under grant agreement n 257462 HYCON2 Network of excellence Chapter 1 Introduction The PnPMPC toolbox is a GNU licensed MatLab toolbox for the modeling of constrained Large Scale Systems LSS with linear time invariant dynamics and for the implementation of the Plug and Play PnP Decentralized De and Distributed Di Model Predictive Control MPC schemes described in 1 2 3 and 4 and for the implementation of large scale estimators and PnP state estimators 5 6 and 4 The PnPMPC toolbox offers also several functionalities for handling zonotopes set and for computing invariant sets The realization of the toolbox has received funding from the European Union Seventh Frame work Programme FP7 2007 2013 under grant agreement n 257462 HYCON2 Network of excel lence The last version of the PnPMPC toolbox can be downloaded at the webpage http sisdin unipv it pnpmpc pnpmpc php Please send bug reports questions or comments to pnpmpc_toolbox unipv it 1 1 Notations e R is the set of real number e N is the set of integers e a bis the
28. addCentralizedInput method The addCentralizedInput method is used to set matrices Been and Deen The method declaration is objlss objlss addCentralizedInput Bcen Dcen i J and it can be used in the following two ways e objlss objlss addCentralizedInput Bcen Dcen in this case we set the global Bcen R and Dcen R matrices of the LSS according to 2 7 and 2 8 e objlss objlss addCentralizedInput Bcen Dcen i j in this case we set Beeni and Deen of matrices Been and Deen 2 5 and 2 6 i e the relation among subsystem i and the centralized input j defined 2 5 and 2 6 The relation among all the other subsystems different from i and centralized control input j is automatically set as zero if not already defined exampleilss 4 specify how exogenous input number 1 affects system i creation of M modelCart modelCart addExoInput 0 1 1 1 4 specify how exogenous input number 1 affects system 2 creation of M modelCart modelCart addExoInput 16 2 2 1 4 specify how exogenous input number 2 affects system 2 creation of M modelCart modelCart addExoInput 0 3 2 2 instead of calling addExoInput three times we can use modelCart modelCart addExoInput 0 0 1 0 16 0 2 3 specify how centralized input number 1 affects system 1 creation of Bcen modelCart modelCart addCentralizedInput 0 4 1 1 specify how centralized input numbe
29. atrices H and K in 2 22 It is also possible to create constraints between two or more subsystems in this case sysindex is a vector of scalars or a cell array of names of the subsystems involved Methods that allow one to add constraints to control inputs outputs exogenous inputs and centralized control inputs objlss objlss addInputConstraint sysindex H K objlss objlss addOutConstraint sysindex H K Note that if the user tries to insert a constraint on the output but the system is in standard setting y x the constraint will be converted in a state constraint and a warning will be displayed objlss objlss addExoConstraint exoindex H K objlss objlss addCentralizedInputConstraint ceninputindex H K 16 and works exactly as addStateConstraint The double methods The double methods return matrices H and K related to constraints The name double has been used for coherency with the MPT2 toolbox 11 The method declaration is Hx Kx objlss doubleStateConstraint index selection and it returns Hx and Kx matrix related to the states of the subsystem subsystems expressed in index index can be a scalar or a string but also a vector of scalars or a cell array of names if one is interested in coupling constraints selection is a Boolean flag with the following meaning 0 exclusive default The state constraints only of the subsystem subsystems indicated in index will be re
30. de the state constraints while fulfilling the matrix inequality A BKaux S A BKaux S lt Q K RKaux aux e S is block diagonal M PCdiag we compute the decentralized symmetric positive definite matrix S and the decentralized static state feedback auxiliary control law Kauxx Kaux diag K1 Km by maximizing the volume of the ellipsoid described by S inside the state constraints while fulfilling the matrix inequality A BKaux S A BKaux S lt Q Kj RKaux e Zero terminal constraint MPCzero we set S 0 and X x 6 2 1 Performance criteria We propose the following performance criteria for evaluating different control schemes e n index Tsim 1 M gt 2 l g k 26 k I where Tsim is the time of the simulation From 6 5 n is a weighted average of the er Qi llug E uf k T airi ror between the real state and the equilibrium state and between the real input and the equilibrium input e index Tetra M 2 5 5 A Prie k Ls 6 6 Taim i 1 JEN where Tim is the time of the simulation and AFie is the power transfer between areas i and j defined in 6 4 This index gives the average power transferred between areas In particular if the 7 index is equal for two regulators the best controller is the one that has the lower value of 39 6 3 Control Experiments We applied the centralized MPC schemes introduced in the previous section to
31. e declare that it is impossible to design a PnPMPC controller for subsystem 7 The interested reader is refer to 1 18 3 2 19 and 4 5 1 1 The pnpmpc method Assume subss stores the model of subsystem 7 given by 2 3 and 2 4 Controller Cj is created through objpnpmpc pnpmpc subss N Xj Uj flag options where N is the receding horizon Xj is a cell matrix Hxj1 K2j1 H2j2 Kxj2 such that parent states j j M are in the set Xj1 xj1 H 1 j1 lt Kan Xj2 xj2 Hj0tj2 lt Kogo where xjl and x72 are states of two parent subsystems of the local subsystem subss Similarly Uj is a cell matrix collecting matrices describing input constraints of parent subsystems options is used to set the sdpsettings for YALMIP see YALMIP manual for details 8 The input argument flag is a structure with the following possible fields e LPdesign true or false If true the design of local PhPMPC controllers is based on the use of LP procedure proposed in 2 3 and Chapter 6 of 4 hence Z is an RCI set and function F will be evaluated online If false the design of local PaPMPC controllers is based on the solution of a nonlinear optimization problem as proposed in 1 and Chapter 5 of 4 in this case parent constraints must be zonotopes Z is an RPI set and function Ki 24 Ty is a linear function Kj aj Ty 27 e ExoBounded indices of exogenous inputs acting on subsystem subss that are distu
32. e estimators All methods can be found in the 1se directory and their utility is described in the help of each function 9 2 Examples In this section we propose some examples Section titles corresponds to m files that can be found in the examples 1se directory 9 2 1 examplellse m Executing file examplellse an example is generated designing the LSS the large scale state estimator and simulating the convergence of the state estimation 52 Chapter 10 Invariant sets In this chapter we present methods to compute invariant sets For more details we refer to the help of each function First of all we give some definitions Definition 8 Robust Positively Invariant RPI set The set X C R is RPI for x t 1 falt w t w t e W CR if a t X gt f a t w t e X Vw t e W Definition 9 Minimal Robust Positively Invariant mRPI set The RPI set X is minimal if every other RPI X verifies X C X The RPI set X 6 is a 6 outer approximation of the minimal RPIX if x E X gt I zx EX andi Bs 0 x x where for 6 gt 0 Bs v x e R z v lt Definition 10 contractive control invariant set The set X C R is a A contractive set for x t 1 f a t u t u t U C R ifVa t X there exist u t U such that z t 1 AX Definition 11 Robust Control Invariant RCI set The set X C R is an RCI set for x t 1 f x t u t w t u t EU CR
33. e frequency deviation Aw and therefore an increase of the power reference APref 6 1 2 Scenario 2 We consider the power network proposed in Scenario 1 and add a fifth area connected as in Figure 6 2 We will simulate Scenario 2 using the load steps specified in Table 6 4 36 Table 6 4 Load of power APr p u for simulation in Scenario 2 APr means a step of required power hence a decrease of the frequency deviation Aw and therefore an increase of the power reference APyre f 6 1 3 Scenario 3 We consider the power network described in Scenario 2 and disconnect the area 4 hence obtaining the areas connected as in Figure 6 3 We will simulate Scenario 3 using load steps specified in Table 6 5 AP op Xj Aw Figure 6 3 Power network system of Scenario 3 37 0 15 0 20 0 15 0 13 0 20 Table 6 5 Load of power APr p u for simulation in Scenario 3 APr means a step of required power hence a decrease of the frequency deviation Aw and therefore an increase of the power reference APre f We can create lss objects for each scenario running makeScenariosPNS For the first sce nario we create 3 files pnsC1 continuous time Iss object pnsD1 discrete time lss object pnsD1ss discrete time system by system lss object Similar files are created for scenario 2 and 3 6 2 Design of the AGC layer for a power network using MPC The goal of the Benchmark is to design the AGC layer for
34. e function returns the H representation of the mRPI Moreover a Polyhedron object of the mRPI is saved in F_epsilon F doubleHK F F_epsilon doublePG If W is a zonotope object the function returns the G representation of the mRPI Moreover a zonotope object of the mRPI is saved in F_epsilonzZ F doublePG F F_epsilonZ 10 1 2 example2epsilon_mRPI m This example shows how to compute an e mRPI set in the case of constraints on the state variables ie r EX 54 Polyhedron eye 2 eye 2 2 2 2 2 epsilon_mRPI A B K W epsilon X 10 2 A contractive control invariant sets In this section we explain how to use the localControlLyapunov class in order to compute a A contractive control invariant set for a discrete time LTI system Our implementation is based on the algorithm proposed in 20 We illustrate the use of the class localControlLyapunov in the following example 10 2 1 examplellocalControlLyapunov m matrices of the LTI system 11 01 J D 0 gd 15 matrices for constraints on state and input variables cx eye 2 eye 2 dx 221 5 2 cu i 3 i J du ES 4 3 dig 4 parameter of the algorithm gt controllability index a parameter of the algorithm we consider the vertices of polytope x cx x lt dx 2 1 5 2 1 5 2 2 2 create the lambda contractive control invariant set localControlLyapunov A B k cx dx cu du x0 plot the lambda contractive control i
35. e given in input using Xj that is a cell matrix Hxj1 K2j1 Hajo Kxj2 such that par ent states zij 7 N are in the set Xj1 xj1 Hxj1xj1 lt Koji Xj2 xj2 Hxj0xj2 lt Kzj2 where xj1 and xj2 are two parent subsystems of the local subsystem subss ctrl is the PnPMPC controller designed with options given in flag 7 1 2 The uk method The control action uj is computed as in 7 1 executing the function u ob jCtrl uRH x0 xj e x0 is the initial local state j 0 a vector of dimension n x 1 e xj is the vector of states from parent subsystems such that K 0 Indices of these parent subsystems are in objCtrl Nij As output we have e u control input value 7 2 Design of unconstrained PnP controllers for a large scale system PnP controllers C for each subsystem in an Iss object are created with the method createCtr1PnP described next 7 2 1 The createCtrIPnP method This method acts on an Iss object and it is called as follows ctrl K createCtrIPnP objlss flag where flag has the same meaning as in the pnpctrl method This method extracts a subsystem calls the constructor of class pnpctr1 and design the controller The output ctrl is an array of dimension 1xnumSys and in the i th position the controller for subsystem is stored K is the overall matrix for the LSS collecting matrices K ij 44 7 3 Examples In this section we propose some examples Section titles corresponds to m files t
36. e three examples that can be found in the PnPMPC toolbox 12 3 1 examplelouterApproximation m definition of function f set X and sampling of f X npoints 500 X zonotope zeros 2 1 2 0 0 3 xf zeros dimension X npoints f x w 1 0 1 x 1 0 5 x 2 exp 0 1 x 1 72 0 1 0 9 x 1 0 1 x 2 0 1 cos x 2 0 05 x 2 72 for i 1 size xf 2 randpoint X xf i f xx XxX end options for outerApproximation function options split ngenerators 2 0 53 max_zono 5 options split alpha options split compute zonotope approximation of f X Zu approximation using 1 zonotope ZuVec approximation using options split max_zono zonotopes 4 Xs zonotope X splitted J jacobian and H hessian are outputs variables The third and fourth input argument are empty not know convex functions gf and hf such that function outerApproximation computes gf and hf for future computations that means the user does f gf hf Therefore the We highlight that the IntLab toolbox http www ti3 tuhh de rump intlab is needed Zu ZuVec Xs J H plots figh figure 1 subplot 1 2 2 hold on Zu plotZ y ZuVec plotZ r plot zf i2 xf 25405 da box on subplot 1 2 1 Xs plotZ g 67 outerApproximation X f options Figure 12 5 Results of example examplelouterApproximation m 12 3 2 example2outerA pproximation m definition of function f
37. er 6 4 Supporting MatLab files In terms of PnPMPC controllers running file makePnpmpcControllersPNS type we can create controllers with different terminal regions and features The input type is a number from 1 to 6 1 Create decentralized PnPMPC controllers based on LP design and equipping each subsystem with input constraints 2 Create decentralized PnPMPC controllers based on LP design and without input constraints 3 Create distributed PnPMPC controllers based on LP design and without input constraints 4 Create decentralized PnPMPC controllers based on nonlinear design and equipping each subsystem with input constraints 40 5 Create decentralized PnPMPC controllers based on nonlinear design and without input constraints 6 Create distributed PnPMPC controllers based on nonlinear design and without input con straints Then we can run each simulation executing file scenario number scenario Full Zero For ex ample if we want run the simulation for Scenario 2 using PnPMPC controllers with zero terminal constraints we should run scenario2Zero options then a file with all simulation data is created We note that using PnPMPC controllers we refer to full when using a quadratic terminal region for controller Cy 1 M Moreover we consider discretization system by system only options is a sdpsettings object for YALMIP The data of a simulation are saved in a file scenario number scenario Full Zero _si
38. er for a LSS The interested reader is refer to 12 13 14 and 15 Details of the design procedures for robust tube MPC can be found in 16 and 17 The model of the LSS is given by 2 7 and 2 8 Our control design procedure will hinge on the following assumptions Assumption 1 The matrix pairs A B Been is controllable Assumption 2 Constraints X and U i M are compact and convex polytopes containing the origin in their nonempty interior 4 1 1 The cenmpc method Assume objlss stores the model of an LSS given by 2 7 and 2 8 MPC controller is created through objcenmpc cenmpc objlss N flag options where N is the receding horizon options is used to set the sdpsettings for YALMIP see YALMIP manual for details 8 The input argument flag is a structure with the following possible fields e ExoBounded indices of exogenous inputs acting on the LSS objlss that are distur bances therefore cenmpc design a robust MPC controller based on tube as in 17 objcenmpc is the MPC controller designed with options given in flag We propose three methods to compute the terminal region and the terminal penalty We refer the reader to help of methods XFqpmax XFqpmaxDec and zeroTerminal 4 1 2 The uRH method In order to compute the control action u we need to solve the MPC optimization problem The control action is computed executing the function u diagnostics objcenmpc uRH x0 din xrefin urefin options
39. erbose 0 one can also use the following instruction ctrliPnpMpcFromLSmodel objlss createCtrlPnPMPC4lsmodel pnsDissLSmodel 4 N kmin kmax m p sdpsettings verbose 0 zero terminal constraint for each pnpmpc controller ctrlPnpMpcZeroFromLSmodel ctrlPnpMpcFromLSmodel for i l size ctrlPnpMpcFromLSmodel 2 ctrlPnpMpcZeroFromLSmodel i ctrlPnpMpcFromLSmodel i zeroTerminal Q R end 60 Chapter 12 Zonotope class In PnPMPC toolbox the zonotope class has been developed as an inherited class of the Polyhedron class of MPT3 7 In the following we describe how to generate zonotope sets Since several func tions have the same meaning of functions of Polyhedron class we defer the reader to the MPT3 manual These functions implement the standard operations between zonotope sets Minkowski sum Pontryagin difference intersection N union U and relational operators C C D D and Zonotope arrays are managed as Polyhedron arrays in MPT3 In the following sections we introduce some useful operators for zonotope sets Most of the definition are from 30 12 1 Creating a zonotope A zonotope in PnPMPC toolbox is created by a call to the zonotope constructor Z zonotope p G This instruction creates a zonotope Z centered in p and described by generators G Example di id 124 3141 zonotope p G plotZ The zonotope is shown in Figure 12 1 In oder to access the G
40. h a zonotope of reduced complexity i e a reduced number of line segment generators Given the zonotope Z the reduction operator redngen Z produces a lower complexity zonotope generated by a maximum of ngen line segment generators The procedure consists of first sorting the columns of G with respect to decreasing Euclidean norm G 91 92 3 Jm Igill 2 lgi all t 1 m 1 Then denoting by Gngen the matrix describing redngen Z we define Gjen G ifm lt ngen Gngen irii Gngen n Gr Gr rS Qngen n 1 9m ifm gt ngen It is important to mention that Z C redngen Z Figure 12 4 shows the application of the reduction operator As one can see a reduction in the number of generators yields a more conservative zonotope The operator redngen has been implemented in PnPMPC toolbox in the following instruction Zo Z reduce ngen The following code shows an example where the number of generators of zonotope Z is reduced Ci 1 1242 45 314561 zonotope p G ngen 3 Zo Z reduce ngen figure 1 hold on Zo plotZ r Z plotZ b 12 2 4 The support function of polytope zonotope sets The support function of a polyhedral set is defined as i 12 1 sup max c g 12 1 One can compute sup using the following instruction sup Z extreme c where X is a polytope or zonotope object 65 20r I I I I I I J 10 5 0 5 10 15 20 15 5 5 5 2 x Figure 12 4
41. hat act on system i and Been ij R and Deen ij R The difference between wy and ucenfi is that up is a local control input i e the output of a local regulator for subsystem while ucen is a global input that cannot be computed locally From models 2 5 and 2 6 the collective dynamics of the resulting LSS is x Ax Bu Md BeenUcen 2 7 y Cx Du Nd Deentcen 2 8 where x zi Tj co X s ER 2 9 is the overall state and n J jem nis u Un Uy e U s e R 2 10 is the overall control input and m em Mi y Un Yo lt gt Ys E R 2 11 is the overall output and p em Pi deR 2 12 collects the exogenous inputs acting on the overall system and Ucen R 2 13 collects the centralized centralized inputs acting on the overall system Moreover one has Au suna Ais B kae Bis A a IJ eR B i e R 2 14 Agi sea Agg Bs Bag All other matrices in 2 7 and 2 8 have a similar block structure Moreover has C R D RP Me R 4 N R 4 Been E R and Deen R The matrix A can also be split into matrices Ag and A defined as An 0 0 Ap 0 6 Ac A Ap 2 15 0 0 Ass Apparently A p collects all the local dynamics of the subsystems while Ac collects coupling terms between subsystems The same decomposition can be done for the matrices B C D appearing in 2 7 and 2 8 Each subsystem can be equipped
42. hat can be found in the examples pnpctrl directory 7 3 1 examplePnpCtrl m We consider a large scale system composed of N masses coupled as in Figure 7 1 vee eve ve 1 1 72 1 TN 1 1 TIN 1 Figure 7 1 Array of masses details of interconnections Each mass M 1 M is a subsystem with input uj and state variables x 24 1 2 i 2 where 21 1 is the displacement with respect to a given equilibrium position x 9 is the horizontal velocity and 100u is an external force in the horizontal direction The values of m 10 while spring and damping coefficients are identical and equal to 0 5 We obtain subsystems Vj by discretizing continuous time models with 0 2 sec sampling time using zero order hold discretization for the local dynamics and treating 7 M as exogenous signals We can successfully design PnP controllers However as the value of masses m i M decreases coupling terms increase and at the same point it becomes impossible to design decentralized controllers Ci For example if all masses m 0 01 we cannot fulfill PnP conditions However by setting dij 1 j N we can completely remove the coupling terms and therefore the synthesis of controllers C amounts to the synthesis of a state feedback controller for each mass without coupling This means that for the system in Figure 7 1 PnP design of distributed controllers Cy is always possible Executing file examplePnpCtrl
43. have the same meaning as in the 1ss class We also highlight that addSystem could also receive one input argument only it can be a subss object see Section 3 3 or a ss object state space object of Control System Toolbox or a tf object transfer function object of Control System Toolbox The addCoupling method By using addSystem method only the block diagonal part of the matrices A B C and D in 2 7 and 2 8 can be set If subsystem 7 is coupled to subsystem j one or more blocks among Ajj Bij Cij or Dij are non zero The method addCoupling allows one to set them The method declaration is objlss objlss addCoupling i j Aij Bij Cij Dij where i and j represent indexes of the the coupled subsystems It is possible to use subsystem names instead of numbers i and j Code for modeling system 3 1 With the three methods presented above it is possible create the model of coupled masses 12 mass elastic constant of the springs 4 damping coefficient 4 creation of an empty lss object modelCart 1ss 1 k M h M B 0 100 M 4 system 1 modelCart modelCart addSystem A B 4 system 2 modelCart modelCart addSystem A B 4 coupling among subsystems 1 and 2 Aij 0 0 K M h M modelCart modelCart addCoupling 1 2 Aij modelCart modelCart addCoupling 2 1 Aij 4 plot graph of subsystems modelCart plot In general the declaration of the methods contains many
44. i th position the controller for subsystem is stored After the creation of the controllers we can run XFqpmax or zeroTerminal methods for setting terminal constraints and then the uRH method to compute the control inputs 5 3 Examples As an example of real application in Chapter 6 we describe PnPMPC controllers for power network systems and run simulations using methods proposed in this chapter In this section we propose simpler example Section titles corresponds to m files that can be found in the examples pnpmpc directory 5 3 1 examplelpnpmpc m We consider a large scale system composed of N masses coupled as in Figure 5 1 ge 1 1 72 1 TN 1 1 TIN 1 Figure 5 1 Array of masses details of interconnections Each mass i M 1 N is a subsystem with state variables z 1 j 2 and input up Ufi 1 Where x 1 is the displacements of mass 7 with respect to a given equilibrium position x 9 is the horizontal velocity of the mass i and 1004 1 is the force applied to mass i in the horizontal direction The values of m have been extracted randomly in the interval 5 10 while spring constants and damping coefficients are identical and equal to 0 5 Subsystems are equipped with the state constraints 21 1 loo lt 1 5 lz gllo lt 0 8 1 M and with the input constraints u j lt Bi where 6 have been randomly chosen in the interval 1 1 5 We obtain models Xj by discretizing co
45. ify norm for the computation of matrices L similar as in 5 and Chapter 3 in 4 ExoBounded indices of exogenous inputs acting on LSS objlss that are distur bances estimator is the large scale state estimator designed with options given in flag 9 1 2 The localEstimator method The state estimation j is computed as in 9 1 executing the function xit yit yitf localEstimator obj i xit xjt yi ui yj d ucen e i is the index of the i th subsystem e xit is Zj at time t e yi is yy at time t e ui is wy at time t e xjt is a vector of x j j E N at time t e yj is a vector of yj 7 E N at time t needed if dij 1 e ucen is Ucen at time t dis d at time t i e exogenous inputs that does not act on the subsystem as not measurable disturbances As output we have e xit is in at time t 1 e yit is yj estimated at time t e yitf is yj estimated at time t 1 9 1 3 The checkErrorEstimation method This method checks if a given initial error e 0 Ei and hence if the state estimator Zi can guarantee convergence of the state estimation and bounded error at all time instants test estimator checkErrorEstimation e i options 5l 9 1 4 The lse2ss method Given a large scale state estimator returns a state space object sys estimator lse2ss i sys is an ss object Note that several methods are available for different features of large scale stat
46. ios are given in Table 6 2 We highlight that all parameter values are within the range of those used in Chapter 12 of 24 We define M as the number of areas in the power network For each scenario discrete time models Xp with Ts 1 sec sampling time are obtained from Di using two alternative discretization schemes e Exact discretization of the overall system acronym D e Discretization system by system i e exact discretization for each area treating ur APrL and aj j Ni as exogenous inputs acronym Dss see 10 for more information In particular we note that Dss preserves the input decoupled structure of DE while D does not 6 1 1 Scenario 1 We consider four areas interconnected as in Figure 6 1 We will simulate Scenario 1 using the load steps specified in Table 6 3 35 EE io fs s w oo 04 10 0 9 0 4 0 1 os Areas mar mas Ilza yllo lt 0 1 zello lt 01 zg alleo 0 1 epa aglloo lt 0 1 l2p 1lloo lt 0 1 lumllo lt 0 5 ugl lt 0 65 lutlo lt 0 65 llull lt 0 55 Ilusio lt 0 5 PR 4 P 2 Pra 2 P 3 Pa5 3 o3 os 1 1 or fo Table 6 2 Model parameters and constraints for systems i 1 5 AP 2 APr 3 APref 2 Awe AEref3 Po Pr P34 APra APLaA Figure 6 1 Power network system of Scenario 1 Table 6 3 Load of power APr p u for simulation in Scenario 1 APr means a step of required power hence a decrease of th
47. iscretization one has to compute the exponential of matrix A Since A is usually a large matrix with many zero elements one has that e computing the exponential is time consuming e elements that are zeros in the matrix A of the continuous time model can be different from zero in the exponential matrix Therefore exact discretization creates coupling between subsystems that were originally decoupled To avoid these problems we also implemented system by system discretization 10 That means that each subsystem in 2 5 and 2 6 is discretized considering states j j 7 i as exogenous inputs We can save computational time because we discretize sequentially subsystems that usu ally have low dimension and we do not create new couplings among subsystems For using this technique set subsystem as method Other properties There are some other properties of 1ss objects that we have not explained yet e couplingmatrix is a matrix with size numSys x numSys composed by boolean elements element i j is true only if subsystem is dynamically coupled or input coupled to sub system j i e i j e name is a cell array of strings with size 1 x numSys that associates a name to every subsystem in the lss object It is possible to set name when calling the addSystem method 10 If a name is not supplied to addSystem a default name will be chosen Every method that needs a subsystem as input has been designed to use these
48. l not be saved sub will have only the couplingA and Aij matrix because the subsystem is dynamically coupled with 236 and 238 so sub couplingA 236 238 sub Aij A237 236 A237 238 where A937 236 is the dynamic coupling among subsystems 237 and 236 and A237 238 is the dynamic coupling among subsystems 237 and 238 3 4 List of other useful methods of class lss e addLocalInput defines a new input for the i th subsystem of an 1ss object e display display properties of an 1ss object in the command window e eig compute eigenvalues of an LSS e generateCouplingmatrix recompute the property couplingmatrix of an lss ob ject get methods return properties of an 1ss object e isCoupled check if two subsystems are coupled e join join two lss objects e name2index index of the subsystem called name e order order of an LSS 22 e plot plot graph of the network of subsystems e plotU plotUcen plotX plotY using the plot function plot signals that represent control inputs centralized control inputs states and outputs e stairsU stairsUcen stairsX stairsY using the plot function plot signals that represent control inputs centralized control inputs states and outputs e set methods set properties of an 1ss object e union join two or more subsystems 23 Chapter 4 The cenmpc class 4 1 Introduction In this chapter we introduce the cenmpc class in order to design an MPC controll
49. le5lss m All the methods described above are useful to add subsystems couplings and constraints to an lss object But if the user inserts wrong parameters it can be useful to remove subsystems features in a selective way This is the goal of the methods discussed next The removeCoupling method objlss objlss removeCoupling i j how This method allows one to remove coupling between subsystems i and j The variable how can be a or b or c or d and specifies if we want to remove coupling terms in the matrices A 18 B C D respectively If how is not passed blocks i 7 in all matrices A B C D will be removed See help removeCoupling for more information The removeExoSignal method objlss objlss removeExoInput i This method allows one to remove the exogenous input dy It deletes the column i of matrices M and N and also the constraints related to this exogenous inputs if they exist Note that a coupling constraint related to exogenous inputs 7 and j after the removal of the i th signal becomes a local constraint of the j th signal For example let du 2dj2 lt 10 be a constraint among exogenous signals dj and dp the execution of objlss objlss removeExoSignal 1 removes dj in the constraint i e the new constraint is 24 lt 10 If one wants to remove all constraints where signal dj appears the removeConstraint method can be used The removeCentralizedInput method objlss objlss remo
50. leCenInputConstraintOnSystem index Next we provide the code for adding constraints described in the code comments to the modelCart object 17 example2lss 0 x_3 1 x_4 lt 3 modelCart modelCart addStateConstraint 2 0 1 3 A 2 x_3 1 x_4 lt 0 and 1 x_3 3 x_4 lt 5 modelCart modelCart addStateConstraint 2 2 1 1 3 0 5 O0 x_1 1 x_2 1 x_3 2 x_4 lt 0 amp 3 x_1 2 x_2 0 x_3 1 x_4 lt 5 modelCart modelCart addStateConstraint 1 2 0 1 1 2 3 2 0 1 0 5 hi 2022 23 modelCart modelCart addInputConstraint 2 2 3 AG 1 F223 modelCart modelCart addInputConstraint 1 2 1 2 3 Hx Kx modelCart doubleStateConstraint 1 disp Hx and Kx are empty because there is no constraint affecting system 1 only Hu Ku modelCart doubleInputConstraint 2 1 disp Hu and Ku contain 2 u_2 lt 3 and i u_1 2 u_2 lt 3 because selection is 1 and we consider all the constraints related to system 2 Input increment In many systems rather than constraining the input it is preferable to bound the input increments dux Upg k 1 upj k With addDeltaUConstraint is possible to define matrices Hsu and K u giving rise to the constraints Hsu d j lt Ksu The method declaration is objlss objlss addDeltaUConstraint sysindex Hdeltau Kdeltau And the corresponding double method is Hdeltau Kdeltau objlss doubleDeltaUConstraint index selection 3 2 5 examp
51. llows one to remove the subsystem number i from an lss object Note that if subsystem 7 is removed this triggers the following changes e if an exogenous signal acts only on subsystem i the method removeExoSignal will be automatically called for removing the signal from the Iss model and a warning appears on the screen The same remark applies to centralized inputs 19 e if two subsystems 7 and j are involved in a coupling constraint and system is removed the constraint becomes a local constraint of j For example let uj up lt 20 be a constraint among subsystems 1 and 2 the execution of objlss objlss removeSystem 1 re moves ujj in the constraint i e the new constraint is up lt 20 In this example we illustrate the use of the removal methods We start from the model created in example4lss example4lss 4 remove exogenous signal number 2 modelCart modelCart removeExoSignal 2 4 remove centralized input number 1 modelCart modelCart removeCentralizedInput 1 remove coupling constraint on the state of system 1 and 2 modelCart modelCart removeConstraint 1 2 state remove dynamic coupling between 1 and 2 so A12 zeros modelCart modelCart removeCoupling 1 2 a remove subsystem number 2 modelCart modelCart removeSystem 2 3 2 6 Performances of the methods dedicated to modeling We discuss here performances in terms of memory occupation and execution time of
52. m type simulation PNS where type simulation depends on the PnP controller design In particular type simulation could be a number from 1 to 4 1 In the simulation decentralized local controllers do not receive state of parent subsystems 2 In the simulation if decentralized local controllers have been designed with LP design we receive state of parent subsystems as in Section 6 5 in 4 3 In the simulation distributed local controllers do not receive state of parent subsystems 4 In the simulation if distributed local controllers have been designed with LP design we receive state of parent subsystems as in Section 6 5 in 4 We highlight that different performances can be achieved using different solvers We also highlight that each simulation could take some minutes to be completely executed If you want take less time for the execution without numerical errors we recommend to install CPLEX optimizer 22 that is free available for academic use We note that most of the functions can be used in a parallel fashion then the interested user can run MatLabpool before the simulations For our simulations results we refer the interested reader to 3 One can also execute all files for modeling and designed PnPMPC controllers running makeAl1PNSfiles type and can delete all files for PnPMPC controllers running clearAl1PNSfiles For each control experiment we provide a file mat of the simulation that contains e lss object
53. method Using zeroTerminal we compute a zero terminal constraint as proposed in Chapter 2 in 15 In particular we consider 64 amp pq k vi A pa k ET E alk ET A olk on k Roya amp vn k Ve pq Ni 0 ary Ni where Q Q gt 0 R and R R gt 0 R Bry k and vag k are the state and input reference trajectories for tracking capabilities We can design the terminal penalt and the ellipsoidal terminal constraint executing the function objpnpmpc objpnpmpc zeroTerminal Q R 28 5 1 4 The uRH method In order to compute the control action uj we need to solve the optimization problem 5 3 The control action is computed executing the function u diagnostics objpnpmpc uRH x0 xj xjinv din xcrefin vrefin options e x0 is the initial local state x j 0 a vector of dimension n x 1 e xj is the vector of states from parent subsystems such that Ki 0 Indices of these parent subsystems are in objpnpmpc Nij e xjinv is a cell array where each element is the state of a parent subsystem Useful only if LPdesign true and there are no disturbances This input will be used to evaluate function Kilt Tuj xy jew for more details see Chapter 6 of 4 e din is the vector of exogenous inputs which act on the subsystem and they are not distur bances Indices of these exogenous inputs are in ob jpnpmpc Exo e xcrefin means the reference state of nominal sy
54. mincon used to solve nonlinear opti mization problems e MPT3 7 which implements the Polyhedron class e YALMIP 8 which include the optimization function solvesdp that is called for solving optimization problems e GraphViz4MatLab toolbox 9 that allows one to plot the graph of a large scale system composed by interconnected systems e INTLAB that allows one to manage interval sets e Symbolic Math Toolbox to manage symbolic variables Please note that GraphViz4MatLab get automatically installed when installing PnPMPC One can check if all required toolboxes are correctly installed with the following instructions yalmiptest mpt_init Also note that part of the toolboxes are needed only for few functions therefore you could not need them 1 3 Installation of the PnPMPC toolbox e Step 1 Add the folder of the PnPMPC toolbox and its subfolders in the MatLab path e Step 2 Run pnpmpc_toolbox_init Remark 1 we tested PnPMPC toolbox on MatLab R2009b and newer version for Windows Linux and MacOS For older versions of MatLab we cannot guarantee that everything works correctly 1 4 Directories The PnPMPC toolbox consists of the following directories e pnpmpc_toolbox cenmpc contains methods for designing MPC controllers e pnpmpc_toolbox epsilon_mRPI contains methods for computing outer approximation of minimal robust positively invariant sets e pnpmpc_toolbox localControlLyapunov contains methods for comp
55. n order to design a state feedback PnPMPC controller for a subsystem PnPMPC controllers are described in 1 18 3 2 19 and 4 In those papers the authors propose a PnP design procedure hinging on the notion of tube MPC 16 and 17 for handling coupling among subsystems and aim at stabilizing the origin of the whole closed loop system while guaranteeing satisfaction of constraints on local inputs and states The model of subsystem i M is given by 2 3 Our control design procedure will hinge on the following assumptions Assumption 3 The matrix pairs Aii Bi Vi M are controllable Assumption 4 Constraints X and U i E M are compact and convex polytopes containing the origin in their nonempty interior Let Z be an RCI set for the nominal subsystem i where the coupling term ew AijXj BijU and bounded disturbances M D are treated as a disturbance We design the following controller for subsystem i Cray Uu Yi Kilti T 5 dig Kig thy 5 1 GEN where K Zi U is any feedback control law guaranteeing Lay Zi gt xt a Zi Va E Xj j E Ni Kij R i and i 0 1 Note that if 6 0 Vi M Vj M the control scheme is completely decentralized since uj depends on local quantities only Following 16 and 17 we set in 5 1 5 2 1Robust Control Invariant RCI set Consider the discrete time Linear Time Invariant LTI system x t 1 Ax t Bu
56. names as an alternative to subsystem numbers However it is better to use numbers when possible because MatLab is slower when working with strings e namex nameu namey named nameucen are cell arrays of strings each with dimen sion 1 x numSys The element in position stores information related to i th subsystem For example namexi is a cell array of n strings which stores the name associated to every state of i th subsystem Similarly remarks can be applied to properties nameu namey named nameucen that are related to local control inputs outputs exogenous inputs and centralized control inputs Names can be written with the ATFXsyntax and are also used in plots of subsystems or time evolution of individual variables 3 2 How to use the lss class through examples In the following sections we provide several examples with comments that illustrate the modeling process with the PnPMPC toolbox Section titles corresponds to m files that can be found in the examples lss directory 3 2 1 examplellss m We show how to model the two coupled mass spring damper systems represented in Figure 3 1 71 1 72 1 Figure 3 1 Array of 2 masses details of interconnections The model is aii T2 wii UT B 3 1 where zp is composed by 1 1 the displacement of first mass with respect to a given equilibrium position and 21 9 the velocity of the first mass Similarly xj is composed by 22 1 and 212 9
57. notopic robust control invariant set Z plot Other useful methods for the zonotopeRCI class are isinside Test if a point is inside the zonotope robust control invariant set Z 1 1 Z isinside x ulnv Given x X compute the control action u x such that x t 1 X Vw W zeros 2 1 Z uInv x 58 Chapter 11 PnPMPC and others toolboxes 11 1 Integration with the toolbox WIDE 29 11 1 1 The lsmodel2lss function In WIDE a large scale system is described by an LSmodel object We use the function lsmodel21ss to convert an LSmodel object to an 1ss object The function declaration is objlss lsmodel21ss objlsmodel m p The arrays m and p are used to select the external inputs and external outputs in the LSmodel object The k th element of the array m is an integer number to classify the k th external input of the LSmodel object if m k 0 the k th external input is a centralized input for the 1ss object if m k 1 the k th external input is an exogenous input for the 1ss object and if m k i the k th external input is a local input for the i th subsystem of the 1ss object Similarly the k th element of the array p is an integer number to classify the k th external output of the LSmodel object if p k i the k th external input is a local input for the i th subsystem of the 1ss object 11 1 2 The createCtrIPnPMPC4lsmodel function Given an LSmodel object we can design PnPMPC controllers u
58. ntinuous time models with 0 2 sec sampling time using zero order hold discretization for the local dynamics and treating j 7 M as exogenous signals At all time steps t the control action uj t computed by the controller C for all i M is kept constant during the sampling interval and applied to the continuous time system Executing file examplelpnpmpc the following functions will be executed In order to create the large scale system composed by N masses we use the function makeModelTrucks1D as follows N 5 modelTrucksiDc modelTrucksiDd makeModelTrucks1D N where N is the number of trucks modelTrucks1Dc and modelTrucks1Dd are the Iss objects of the continuous time and discrete time models respectively in form of lss objects We can design local controllers executing the function 30 options sdpsettings verbose 0 TrucksiDpnp makePnpmpcTrucks1D modelTrucksiDd options whose instructions are described next More in detail the controllers Cj i M are generated as follows Trucks1Dpnp is an array of N pnpmpc objects NMPC 10 prediction horizon TrucksiDpnp modelTrucksiDd createCtrlPnPMPC NMPC Then we add a zero terminal constraint to each controller using the following instructions 10 eye 2 R 1 parfor i 1 modelTrucksiDd numSys TrucksiDpnp i Trucks1iDpnp i zeroTerminal Q R end We note that the functions makePnpmpcTrucks1D or modelTrucks1Dd createC
59. nvariant set plot Contractiveness lambda L lambda Other useful methods for the localControlLyapunov class are isinside Test if a point is inside the A contractive control invariant set L 1 51 L isinside x 55 ulnv Given x X compute the control action u x such that x t 1 AX zeros 2 1 L ulInv x double The function returns the H representation of a A contractive control invariant set More over a Polyhedron object of the set can be obtained through the attribute LP L double LLP 10 3 Parameterized robust control invariant sets In this section we explain how to use the parameterizedRCI class in order to compute a parameterized robust control invariant set for a linear constrained systems Our implementation is based on the algorithm proposed in 20 We illustrate the use of the class parameterizedRCI in the following example 10 3 1 examplelparameterizedRCI m matrices of the LTI system 11 01 1 C 0 11 matrices for constraints on state and input variables eye 2 eye 2 1 2 2a 2 Ci 1 Ty 3 3 4 parameter of the algorithm greater than the controllability index of A k k 5 x0 is a parameter of the algorithm such that W subseteq hull x0 x0 0 99 0 99 0 99 0 01 0 01 0 01 0 01 0 99 create the parameterized robust control invariant set parameterizedRCI A B k cx dx cu du x0 plot the parameterized robust
60. omputer Aided Control Systems Design Taipei Taiwan September 2 4 2004 pp 284 289 M Dunham K Murphy L Peshkin and D Eaton GraphViz4Matlab 2004 Online Available https github com graphviz4matlab graphviz4matlab M Farina P Colaneri and R Scattolini Block wise discretization accounting for structural constraints Automatica vol 49 no 11 pp 3411 3417 2013 M Kvasnica P Grieder and M Baoti Multi Parametric Toolbox MPT 2004 Online Available http control ee ethz ch mpt 2 A Bemporad and M Morari Control of systems integrating logic dynamics and con straints Automatica vol 35 no 3 pp 407 428 1999 72 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 D Q Mayne J B Rawlings C V Rao and P O M Scokaert Constrained model predictive control Stability and optimality Automatica vol 36 no 6 pp 789 814 2000 J M Maciejowski Predictive Control with constraints Upper Saddle River NJ USA Prentice Hall 2002 J B Rawlings and D Q Mayne Model Predictive Control Theory and Design Madison WI USA Nob Hill Pub 2009 D Q Mayne M M Seron and S V Rakovi Robust model predictive control of con strained linear systems with bounded disturbances Automatica vol 41 no 2 pp 219 224 2005 S V Rakovi and D Q Mayne A simple tube controller
61. optimal parameters that are op tional and the user can provide only meaningful quantity In the example above parameters Ci Di name Ts are not supplied to addSystem method because we are assuming y x so we are in the standard setting discussed in Section 3 1 Moreover we do not want to associate a name to the model and we are defining a model in continuous time so that Ts is empty Similar remarks apply to all methods of the toolbox 3 2 2 example2lss m Now we will add an exogenous input and a centralized control input to the lss object modelCart through examplellss m With reference to 2 7 and 2 8 we want to set op de e di ood e Like in the previous section we start with the description of the relevant methods The addExoInput method The addExoInput method is used to set matrices M and N The method declaration is objlss objlss addExoInput M N i 4 13 and it can be used in the following two ways e objlss objlss addExoInput M N in this case we set the global M R 4 and N RP 4 matrices of the LSS according to 2 7 and 2 8 e objlss objlss addExoInput M N i j in this case we set only blocks Mi and Nij of matrices M and N 2 5 and 2 6 i e the relation among subsystem i and the exogenous signal j defined 2 5 and 2 6 The relation among all the other subsystems different from i and exogenous signal j is automatically set as zero if not defined yet The
62. r bances therefore pnpmpc design a robust PnPMPC controller as in Chapter 7 of 4 e distributedParents indices of subsystems such that 6 1 in 5 1 It must be a subset of the indices of parent subsystems e norm specify norm for the computation of matrices K see Chapter 7 in 4 e boundKij specify bound for the computation of matrices K see Chapter 7 in 4 e k if LPdesign true k kmin kmax allows to try different k for the LP design procedure proposed in 2 3 and Chapter 6 in 4 e Qi and Ri if LPdesign false Qi and Ri are inputs for the nonlinear optimization problem as proposed in Chapter 5 in 4 objpnpmpc is the PnPMPC controller designed with options given in flag We propose two methods to compute the terminal region and the terminal penalty 5 1 2 The XFqpmax method Using XFqpmax we compute a quadratic terminal region using procedures proposed in 21 In particular we consider Lil k vq k lk kQ llk ET 1 walk un k Riog lk v 6 Vr a Ni a Na E NDV S l t Ni DND 8 R SB lt 1 where Q Q gt 0 R R R gt 0 R and S S gt 0 eR Sre k and vie k are the state and input reference trajectories for tracking capabilities We can design the terminal penalty S and the ellipsoidal terminal constraint executing the function objpnpmpc objpnpmpc XFqpmax Q R options 5 1 3 The zeroTerminal
63. r 1 affects system 2 creation of Bcen modelCart modelCart addCentralizedInput 0 5 2 1 instead of the previous calls to addCentralizedInput we can use modelCart modelCart addCentralizedInput 0 4 0 5 14 plot all signals modelCart plot all If like in this case the user does not need to set the N or Dcen matrix it is possible to use empty matrices Empty or not passed parameters are given default values The same remark applies to all other methods in the toolbox 3 2 3 example3lss m This example shows how to exploit modularity in the subsystem interconnection for quickly cre ating an LSS The code below allows one to create a string of N mass spring damper systems as shown in Figure 3 2 with different M mass different k elastic constant of the springs and different h damping coefficient verifying the bounds minM lt M lt maxM mink lt k lt maxk and minh lt h lt maxh U2 U N 1 gt eve 1 1 72 1 TN 1 1 TIN 1 Figure 3 2 Large scale system to model in example3lss The system is characterized by the matrices 0 1 0 1 An k 1 h 1 gt ANN k N 1 ufo MA MI MN MN Vi 2 N 1 0 1 0 A ka po _ ha pa gt Bi 100 3 5 M i M i M 0 0 0 0 Ai k i 1 h i 1 gt Ai i 1 k i 1 ne MG l M i 1 Mi Mi Qu 1 o Dy 0 The MatLab code when N 1000 is reported in the follo
64. representation of zonotope Z one can access to properties p and G In oder to access the H representation of zonotope Z one can use the command computeHRep Zo Z computeHRep 61 Figure 12 1 Zonotope Z H Zo A K Zo b Since the computation of the H representation of zonotope Z requires the computation of the V representation both representation are stored in the output variable Zo for future uses 12 2 Functions for zonotope objects 12 2 1 Bisection and Split operators The operator Bisect generates two sub zonotopes from one zonotope In particular given a zonotope Z the operator Bisect Z generates the two sub zonotopes Jk Ik Jk Jk Z p 9 Ola Fo Ild ldlliny lt 1 where gx is the k th column of G and is the biggest generator i e we split the k th column in the middle With the operator Splitk Z a it is possible to split the k th column of G in a desired position a i e Split Z a generates two sub zonotopes Z p gk 1 0 91 gna 9m d ld ling lt 1 Z p gra g1 g9k 1 a 9ml d dllinf lt 1 where gx is a column of G matrix and the parameter a 0 1 Figure 12 2 shows an example of the operator Bisect Z applied to the zonotope Z in Figure 12 1 This example shows that the operator Bisectx Z can generate two sub zonotopes which intersect in their interior The reason of the overlapping of Z and Z is that the line segment generators g1
65. respectively 100uj 9 is the force applied to mass in the horizontal respectively vertical direction The values of m have been extracted randomly in the interval 8 10 while spring constants and damping coefficients are identical and equal to 0 5 Subsystems are equipped with the state constraints x loo lt 1 5 31 Figure 5 2 Array of masses details of interconnections j 1 3 zugllo lt 0 8 M 1 2 4 and with the input constraints ujjl o lt Bi where Bi have been randomly chosen in the interval 1 1 5 We obtain models X by discretizing continuous time models with 0 2 sec sampling time using zero order hold discretization for the local dynamics and treating xj j Ni as exogenous signals 10 At all time steps t the control action u t computed by the controller C for all i M is kept constant during the sampling interval and applied to the continuous time system 50 40 Sf 307 20 Figure 5 3 Position of the N 1024 trucks on the plane In order to create the large scale system composed by N masses we can use the function 32 makeModelTrucks2D as follows 16 modelTrucks2Dc modelTrucks2Dd makeModelTrucks2D N where N is the number of trucks modelTrucks2Dc and modelTrucks2Dd are the continuous time and discrete time large scale models respectively in form of 1ss objectes We can design local controllers executing the function T
66. rucks2Dpnp makePnpmpcTrucks2D modelTrucks2Dd whose instructions are described next Controllers Cj M are generated as follows Trucks2Dpnp is an array of N pnpmpc objects NMPC 10 prediction horizon Trucks2Dpnp modelTrucks2Dd createCtrlPnPMPC NMPC Then we add a zero terminal constraint to each controller using the following instructions Q 10 eye 4 R eye 2 parfor i 1 modelTrucks2Dd numSys Trucks2Dpnp i Trucks2Dpnp i zeroTerminal Q R end We note that the functions makePnpmpcTrucks2D or modelTrucks2Dd createCtr1PnPMPC could take some minutes to be completely executed If you want take less time for the execution we recommend to install CPLEX optimizer 22 that is free available for academic use To start a simulation we can use the function runSimTrucks2D in the following way repmat 1 0 1 0 N 1 30 runSimTrucks2D Tsim x0 modelTrucks2Dc modelTrucks2Dd Trucks2Dpnp This function computes the control action using the method uRH of the pnpmpc class We did not specify a solver in order to be as general as possible We really recommend the use of CPLEX In 2 3 and Chapter 6 of 4 we propose the case of N 1024 trucks and we give the results in terms of computational times 33 Chapter 6 Hycon2 Benchmark Power Network System Note for MatLab simulations in this chapter we recommend the use of CPLEX optimizer 22 6 1 Introduction An example of a real
67. scenarios 1 2 and 3 Furthermore for each scenario we discretized the continuos system with both discretization schemes D and Dss At time t we solve the optimization problem 6 3 and then apply the control action to the continuos time system keeping the value constant between time t and t 1 If at time t the power load increases or decreases we assume the controller can use this information at time t This means at time t the controller knows exactly the value of AP hence can use it We highlight that violation of this assumption can impact considerably on the index 7 In all experiments we use Tsim 100 In Table 6 6 and 6 7 the values of the performance parameters 7 and respectively are reported for each control experiment Scenario 1 Scenario 2 Scenario 3 D Da D De D De RPC Full 00219 00249 0056 0 0347 00510 0 051 MPC diag 0 0249 0 0249 0 0346 0 0547 00510 0 051 MPC zero 00259 0 0249 00346 00347 00510 0 051 Table 6 6 Values of the performance parameter n using different centralized MPC schemes for the AGC layer CT Seen Scenario Scenario 3__ p ps p ps 0 Ds PC 00030 0 0029 00063 0 0060 0 0060 0 0055 MPC diag 0 0080 0 0029 00068 0 0061 0 0060 0 0058 MPC zero 00030 0 0025 0 0065 0 0059 0 059 0 0055 Table 6 7 Values of the performance parameter using different centralized MPC schemes for the AGC lay
68. set of integer a a 1 b e A gt 0 means that the matrix A is positive semi definite A gt 0 means that matrix A is positive definite e A diag Aii Ass is the block diagonal matrix An 0 0 0 0 1 1 0 0 Ags Tij Ul Yq di U cen are column vectors of suitable dimensions e x i s means that x is composed by column vectors aj i 1 s stacked in a single column i e tl XQ x 1 2 Ts e The symbol denotes the Minkowski sum e g A B C if and only if A a a b c Vb e B Vc Ch e Di G G19 0G e x indicates the cartesian product and _ Gi G1 x x Gs Definition 1 A polyhedron X is the intersection of a finite number of closed half spaces Therefore we can represent a polyhedron either as X xeR Hr lt K where H R and K R H representation or as a convex hull of vertices V representation Definition 2 A zonotope is a centrally symmetric convex polytopes given a vector pe R and a matrix R the zonotope X C R is the set X x x p Ed d xo lt 1 with d e R We will refer to this representation as G representation 1 2 Required toolboxes For a complete use PnPMPC toolbox requires the following toolboxes to be installed e Control System Toolbox which implements the class ss state space We also use the function c2d for time discretization e Optimization Toolbox which implements the class f
69. sing the createCtr1PnPMPC41smodel function This function calls the 1smode121ss function and the createCtr1PnPMPC method of the 1ss object of the function and returns an array of pnpmpc objects The function declaration ne n ctrl objlss createCtrlPnPMPC4lsmodel objlsmodel N k m p options For the meaning of the input arguments we defer the reader to Section 11 1 1 and Section 5 2 1 11 1 3 Example We create an LSmodel object of the power network system proposed in Scenario 1 in Chapter 6 We can create the LSmodel object executing the file makePNSLSmodel m Executing the file makePnpmpcControllersPNS4LSmodel we can design PnPMPC controllers for the power network described by an LSmodel We show below the essential part of the code diag 500 0 01 0 01 10 10 15 kmin 88 8 6 kmax 20 20 20 20 m and p array for lsmodel21ss the odd elements of m are the power references of each area therefore they are the local inputs of each subsystem the even elements of m are the power loads of each area therefore they are the exogenous signals of each subsystem 11 12 13 1 4 1 the elements of p are the outputs of each subsystem therefore they are the states of each area 1111222233334444 convert LSmodel to lss objlss lsmodel21ss pnsDissLSmodel m p design pnp controllers for i 1 numel kmin flag i k kmin i kmax i end ctrlPnpMpcFromLSmodel objlss createCtr1lPnPMPC N flag sdpsettings v
70. stem e over the control horizon it can be expressed as a vector n x 1 reference constant over the control horizon or a matrix n x N the columns 1 refer to the values at the time instant 1 within the control horizon columns 2 refer to the values at time 2 etc If this parameter is empty or not passed it is set to zero by default e vrefinin means the reference input of nominal system vre over the control horizon it can be expressed as a vector m x 1 reference constant over the control horizon or a matrix mi x N the columns 1 refer to the values at the time instant 1 within the control horizon columns 2 refer to the values at time 2 etc If this parameter is empty or not passed it is set to zero by default e options sdpsettings object for YALMIP options As output we have e u control input value e diagnostics information about the optimization problems 5 2 Design of PnPMPC controllers for a large scale system PnPMPC controllers Cj for each subsystem in an lss object are created with the method createCtr1lPnPMPC described next 5 2 1 The createCtrIPnPMPC method This method acts on an lss object and it is called as follows ctrl createCtrlPnPMPC 0bjlss N flag option 29 where N flag options have the same meaning as in the pnpmpc method This method extracts a subsystem calls the constructor of class pnpmpc and design the controller The output ct r1 is an array of dimension 1xnumSys and in the
71. t w t with x t R u t R w t W and subject to constraints u t U C R and w t W C R The set X C R is an RCI set with respect to w t W if Ve t X then there exist u t U such that x t 1 X Vw t W Model of subsystem i without coupling terms 26 where vp Olt and w 0t are optimal values of the variables vjj 0 and u 0 respectively appearing in the MPC problem Ni 1 Prat min I G vg Vi ga 5 3a Via 0 N 1 k 0 a 0 zalt 14 0 Zi 5 3b p k 1 Aiia k e Bivjg k ke0 Nj 1 5 3c Bray k Xi ke0 N 1 5 3d va k E Vi ke0 Nj 1 5 3e a Ni e X 5 3f In 5 3 N N is the control horizon 4 R R is the stage cost Vy R R is the final cost and Xp is the terminal set Moreover from 5 3c and 5 3e the tightened constraints and V are defined respectively as Therefore in order to design a PnPMPC controller for subsystem i we need to solve the following problem Problem P Given 6 compute matrices K and nonempty RCI Z for the nominal subsys tem 7 treating the coupling term as a disturbance Compute sets X and Vj We highlight that Problem P can be solved using efficient procedure proposed in 20 that requires the solution of a suitable LP problem Moreover X and V can be computed using optimizers from the LP problem If Problem P cannot be solved w
72. tates and disturbances Moreover the subsystems are input decoupled We define for M the Local State Estimator LSE as Sa Fy At Baug L y Cit DI At DL dig li lun Gim JEN JEN gt Bi Uj BeenUcen Mid JEN 9 1 where R is the state estimate Li R Pi are gain matrices and d 0 1 This implies that Zi depends only on local variables jj uj and yj and parents variables ujj Yup IE N centralized inputs cen and exogenous inputs d Binary parameters big GEN can be chosen equal to one for exploiting the knowledge of parents outputs or equal to zero for reducing the number of transmitted output samples The lse class design state estimators Za i M guaranteeing converge of state estimation and bounded error at all time instants The complete algorithm design can be found in Chapter 3 in 4 9 1 1 The lse method Assume subss stores the model of subsystem i given by 2 5 and 2 6 The collection of state estimators Zi is created through estimator lse objlss flag options where e objlss is the large scale system e options is the sdpsettings structure for solvesdp options in YALMIP 50 e The input argument flag is a structure with the following possible fields Delta square matrix of dimension objlss numSys Element ij is false if sub system j does not send outputs to subsystem i i e Lj 0 true otherwise default true norm spec
73. the methods of the PnPMPC toolbox dedicated to modeling As a first case study we use the string of N masses described in examplellss m and example3lss m with local constraints on states and inputs In the left panel of Figure 3 3 one can see e the execution time for creating the lss object blue line e the execution time for discretizing the lss object using the system by system approach red line as a function of the number of masses In the right panel of Figure 3 3 it is reported e the memory occupation in Kb due to the continuous time model blue line e the memory occupation in Kb due to the discretized model red line as a function of the number of masses These graphics have been obtained on a processor Intel Core i3 3 06 GHz with 4GB of Ram 1 33 GHz running MatLab r2011a 20 Execution time Memory occupation 45 T T T T 8000 T T T 7000 6000 5000 4000 t h Mem Kb 3000 2000 1000 Figure 3 3 Performance of PhPMPC toolbox 3 3 The subss class This class allows one to model individual subsystems and it is useful for the following reasons 1 Creating a subss object and pass it to an lss object with the addSystem method In this way dynamics and constraints of the subsystem will be added to the LSS 2 Extracting a subsystem from an lss object through the get System method Many methods of this class have the same name of methods of the lss class and work in a
74. to 6 and Chapter 8 in 4 8 1 1 Bounded error state estimation In several applications we need to guarantee that the state estimation error is bounded i e u u E Methods proposed in 6 and Chapter 8 in 4 allow one to guarantee prescribed bounds 46 In the following we introduce methods needed to equip an 1ss object with state estimation error constraints Methods for adding state estimation error constraints The addErrorEstimationConstraint method sets the matrices He and K such that E Xj fp E R He cu La lt Ke The method declaration is objlss objlss addStateConstraint sysindex H K Where sysindex is the number or name of the subsystem to which constraints have to be added and H and K are the matrices He and K It is also possible to create constraints between two or more subsystems in this case sysindex is a vector of scalars or a cell array of names of the subsystems involved In a similar way as in Section 3 2 4 method He Ke objlss doubleErrorEstimationConstraint index selection E returns matrices associated to sets E 8 1 2 The pnpEstimator method Assume subss stores the model of subsystem i given by 2 5 and 2 6 State estimator Zi is created through estimator pnpEstimator subss Ej Cj 0j flag where e Ej is a cell matrix Hej1 Kej1 Hej2 Kej2 such that parent state estimation errors
75. tomatica vol 44 no 1 pp 216 224 2008 D M Raimondo S Riverso C N Jones and M Morari A Robust Explicit Nonlinear MPC Controller with Input To State Stability Guarantees in Proceedings of the 18th IFAC World Congress Milano Italy August 28 September 2 Aug 2011 pp 9284 9289 74
76. tr1PnPMPC could take some minutes to be completely executed if N is large If you want take less time for the execution we recommend to install CPLEX optimizer 22 that is free available for academic use To start a simulation we use the function runSimTrucks1D in the following way repmat 1 0 N 1 30 runSimTrucks1D Tsim x0 modelTrucks1Dc modelTrucksiDd TrucksiDpnp This function computes the control action using the method uRH of the pnpmpc class We really recommend the use of CPLEX 5 3 2 example2pnpmpc m In this example differently from examplelpnpmpc m we consider that each subsystem is affected by bounded disturbances Therefore in order to design robust PnPMPC controllers we select the bounded exogenous inputs using the struct flag and its field ExoBounded 5 3 3 example3pnpmpc m We consider a large scale system composed by N masses coupled as in Figure 5 3 for the case of N 1024 where the four edges connected to a point correspond to springs and dampers arranged as in Figure 5 2 Hereafter we assume that N z for some z N z gt 1 Each mass i M 1 N is a subsystem with state variables x 4 1 2 2 2 i 3 Tji 4 and input up 1 2 where y and y are the displacements of mass 7 with respect to a given equilibrium position on the plane equilibria lie on a regular grid x 2 and 4 are the horizontal and vertical velocity of the mass i respectively and 100u 1
77. turned 1 all All state constraints which includes the subsystem subsystems indicated in index will be returned For example let c be a constraint related to subsystem 1 c3 be a constraint related to subsystem 3 and c13 be a coupling constraint between subsystem 1 and 3 If index 3 and selection 0 or not passed H and K contain only c3 because this is the only constraint that involves subsystem 3 only If selection 1 H and K will contain c3 and c13 i e all constraints involving subsystem number 3 Another possible use of this method is X objlss doubleStateConstraint index selection with only one output argument In this case X is the Polyhedron i e an object of Polyhedron class see the documentation of the MPT3 for help of subsystem specified in index giving rise to a state constraint Also if index contains more than one index a Polyhedron object will be returned Other double methods with similar meaning and use are Hu Ku objlss doubleInputConstraint index selection Hy Ky objlss doubleOutConstraint index selection Hd Kd objlss doubleExoConstraint index selection Hcen Kcen objlss doubleCenInputConstraint index selection Moreover the following methods return matrices Hd Kd and Hcen Kcen for exogenous inputs and centralized control inputs acting on subsystem index Hd Kd objlss doubleExoConstraintOnSystem index Hcen Kcen objlss doub
78. ulations as scenario number scenario Full Zero _sim type simulation PNS Assume we want to simulate Scenario 2 where decentralized controllers have been designed using non linear design with zero terminal constraint and input constraints We need to load MatLab file scenario2ZeroData_sim1PNS This operation can be performed with the MatLab command load scenario2ZeroData_simiPNS e Step 3 Start a simulation from Simulink will produce the results of the control experiments sim simulatorPNS_AGC_2 Simulation results are summarized in 3 1 25 and Chapters 5 6 and 7 in 4 42 Chapter 7 The pnpctrl class 7 1 Introduction In this chapter we introduce the pnpctrl class in order to design an unconstrained PnP controller for a single subsystem Unconstrained PnP controllers are described in 19 and Chapter 4 in 4 In those papers the authors propose a PnP design procedure hinging on the notion of small gain theorem for networks 26 for handling coupling among subsystems and aim at stabilizing the origin of the whole closed loop system The model of subsystem i i M is given by 2 3 and 2 4 We consider that subsystems are input decoupled Our control design procedure will hinge on the following assumption Assumption 5 The matrix pairs Aii Bi Vi M are controllable We design the following controller for subsystem i JEN where K R s and d 0 1 Note that if d 0 Vi M Vj
79. uting control invariant sets e pnpmpc_toolbox lse contains methods for designing large scale state estimators e pnpmpc_toolbox lss contains methods for modeling of LSS e pnpmpc_toolbox parameterizedRCI contains methods for computing robust con trol invariant sets e pnpmpc_toolbox pnpctr contains methods for designing of unconstrained PnP con trollers e pnpmpc_toolbox pnpEstimator contains methods for designing of PnP state esti mators pnpmpc_toolbox pnpmpc contains methods for designing of PaPMPC controllers pnpmpc_toolbox subss contains methods for the modeling of a subsystem of an LSS pnpmpc_toolbox zonotope contains methods for using zonotope sets pnpmpc_toolbox zonotopeRCI contains methods for computing zonotopic robust control invariant sets pnpmpc_toolbox examples contains examples demonstrating the functionalities of PaPMPC toolbox e pnpmpc_toolbox extra contains extra useful functions e pnpmpc_toolbox toolbox contains additional toolboxes provided with PnPMPC tool box For familiarizing quickly with PnPMPC toolbox examples are supplied in the directory pnpmpc_toolbox examples These m files can also be used as templates for your personal experiments Chapter 2 Modeling of interconnected systems In this chapter we describe the class of systems considered in the PnPMPC toolbox Features of subsystems and their interconnections will be discussed in details The use
80. veCentralizedInput i This method allows one to deletes the i th centralized control input It eliminates the column i of matrices Bcen and Dcen and also the constraints related to this centralized control input if they exist Note that a coupling constraint related to centralized control inputs i and j after the removal of becomes a local constraint of j If one wants to remove all constraints where centralized input cen ij appears the removeConstraint method can be used The removeConstraint method objlss objlss removeConstraint index flag This method allows one to remove constraints index can be as usually the subsystem number or name if we want to remove local constraints But it can also be a vector of scalars or a cell array of names if we want to remove coupling constraints flag can be state or in or out or exo or cen or deltau In this way it is possible to remove for instance state constraints only If flag is exo or cen we refer to the constraints of exogenous input dj or centralized control input cen i and in this case we have to specify the number of signal signals interested If flag is not specified all constraints on states inputs and outputs of the system in index will be removed Fore example if flag is empty and index is 1 then X U and Y will be removed See help removeConstraint for more information The removeSystem method objlss objlss removeSystem i This method a
81. wing script N 1000 minM 1 maxM 10 mink 0 1 maxk 0 9 minh 0 1 maxh 10 M minM rand 1 N maxM minM 15 w Il mink rand 1 N 1 maxk mink minh rand 1 N 1 maxh minh p Il modelCart 1ss modelCart modelCart addSystem 0 1 k 1 M 1 h 1 M 1 0 100 M 1 1 0 0 for i 2 N 1 modelCart modelCart addSystem 0 1 k i 1 k i M i h i 1 h i M i 0 100 MC i 1 0 0 modelCart modelCart addCoupling i 1 i1 0 0 k i 1 M i 1 h i 1 M i 1 zeros 2 1 modelCart modelCart addCoupling i i 1 0 0 k i 1 M i h i 1 M i zeros 2 1 end modelCart modelCart addSystem 0 1 k N 1 M N h N 1 M N 0 100 M N 1 0 0 modelCart modelCart addCoupling N N 1 0 0 k N 1 M N h N 1 M N zeros 2 1 modelCart modelCart addCoupling N 1 N 0 0 k N 1 M N 1 h N 1 M N 1 zeros 2 1 3 2 4 example4lss m This example shows how to add state and control input constraints to an lss model The addition of constraints on output exogenous input and centralized control input can be done in a similar way Relevant methods are described next Methods for adding constraints The addStateConstraint method sets the matrices H and K in 2 22 The method declaration is objlss objlss addStateConstraint sysindex H K Where sysindex is the number or name of the subsystem to which constraints have to be added and H and K are the m

PnPMPC Toolbox v. 1.0 - User manual

Contents

Download Pdf Manuals

Related Search

Related Contents