Home

PnPMPC Toolbox v. 0.9 - User manual

1. 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 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 Xij C RT 2 172 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 19 ieM i jeMi j Moreover we can define constraints for the exogenous signals and for the centralized inputs deDCR 2 20 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. For the meaning of constants as well as some typical parameter values we refer the reader to Table 5 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 percent change in load Defined as 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 5 1 Variables of a generation area with typical value ranges 15 p u stands for per unit We note that model 5 1 is input decoupled since both AP and APr act only on subsystem DE Moreover subsystems DE are parameter dependent since the local dynamics depends on the yen 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 A0 an
3. 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 and a warning appears Every method that needs a subsystem as input has been designed to use these 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 namex 1 i is a matrix of strings of n rows which stores the name associ ated to every state of i th subsystem Similar remarks can be applied to properties nameu namey named nameucen that are related to system input output exogenous signals and centralized input Names can be written with the BTfXsyntax 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 Fig
4. 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 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 lss and subss classes 3 1 The Iss class There are many MatLab toolboxes that offer facilities for modeling MIMO LTI systems However most of them have not been designed to handle the interconnection of several subsystems There fore 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 properties of an lss 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 ou
5. 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 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 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 Automatica vol 44 no 1 pp 216 224 2008 62
6. kmax 20 20 20 20 m and p array for lsmodel21lss 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 Ea 1 2 1 3 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 ctrlPnpMpcFromLSmodel objlss createCtr1PnPMPC N kmin kmax sdpsettings verbose 0 one can also use the following instruction ctrlPnpMpcFromLSmodel objlss createCtrlPnPMPC4lsmodel pnsDissLSmodel 4 N kmin kmax m p sdpsettings verbose 0 4 zero terminal constraint for each pnpmpc controller ctrlPnpMpcZeroFromLSmodel ctrlPnpMpcFromLSmodel for i l size ctrlPnpMpcFromLSmodel 2 ctrlPnpMpcZeroFromLSmodel i ctrlPnpMpcFromLSmodel i zeroTerminal Q R end 49 Chapter 8 Zonotope and Polytope classes In PnPMPC toolbox the zonotope class has been developed as an integration of MPT2 3 New methods have been proposed also for the polytope class In the following we describe how to generate zonotope sets Since several functions have the same meaning of functions of polytope class we defer the reader to the MPT2 manual These functions implement the standard op erations betwee
7. objlss objlss removeCentralizedInput i This method allows one to deletes the centralized input number i It eliminates the column i of matrices Bcen and Dcen and also the constraints related to this centralized input if present Note that a coupling constraint related to centralized inputs and j after the removal of i becomes a local constraint of j If one wants to remove all constraints where centralized input Ucen i 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 to exogenous signals dj or centralized input 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 allows one to remove
8. rs G is a n x n diagonal matrix rs G ii 5 Gi l i Liscia j 1 52 The operator rs has been implemented in PnPMPC toolbox in the following instruction BB bounding_box Z 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 C1 1 13 1242 45 31456 245 1 5 Z zonotope p G BB bounding_box Z figure 1 hold on plot BB plot Zz 7b Figure 8 3 An example of bounding box of a set Z 8 2 3 Reduction operator Another important operator is the reduction operator whose purpose is to outer bound a given zonotope with 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 9ml gil gt Ilgi4all i Long 1 53 Then denoting by Gngen the matrix describing redngen Z we define Ggen G ifm lt ngen Groen 91 lt 3 Ingen n gr G rs Ingen n 1 arr Gm ifm gt ngen It is important to mention that Z C redngen Z Figure 8 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 im
9. 2 8 e objlss objlss addExoInput M N i j in this case we set only blocks M i 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 already defined The 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 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 inthis 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 sub systems different from i and centralized input j is automatically set as zero if not already defined exampleilss 4 specify how exogenous input number 1 affects system 1 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 modelC
10. are in bold 39 Area 1 3 Area2 9 10 4 10 ohj 2 xv 3 3 lt 4 2 0 I 4 2 0 50 100 0 50 100 t s t s 3 Area3 3 Area4 4x 10 5x 10 2 0 e So lt 0 aq _2 l 5 4 10 0 50 100 0 50 100 t s t s a Frequency deviation in each area controlled by the DePnPMPC bold line and centralized MPC dashed line Area 1 Area 2 0 2 0 2 0 15 0 1 57 5 o 0 1 a 0 lt lt 0 05 0 1 0 2 50 100 o 0 50 100 t s t s Area 3 Area 4 0 2 0 4 Pr 4 0 2 z O o 0 o lt lt 0 1 o 0 25 50 100 0 25 50 100 t s t s b Load reference set point in each area controlled by the DePnPMPC bold line and centralized MPC dashed line Figure 5 4 Simulation Scenario 1 5 4 a Frequency deviation and 5 4 b Load reference in each area 40 Area 1 gt 2 Area 2 gt 3 0 01 0 01 0 005 a opf 3 A 2 Ko a o 0 0 01 7 0 005 0 02 0 01 s 0 50 100 a0 0 50 100 t s t s Area 3 gt 4 0 04 sz 0 02 2 ni lt 0 i 0 02 0 50 100 t s Figure 5 5 Simulation Scenario 1 tie line power between each area controlled by the DePnPMPC bold line and centralized MPC dashed line 41 Chapter 6 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
11. frequency deviation Aw and therefore an increase of the power reference APref 5 1 2 Scenario 2 We consider the power network proposed in Scenario 1 and add a fifth area connected as in Figure 5 2 We will simulate Scenario 2 using the load steps specified in Table 5 4 33 Table 5 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 APre f 5 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 5 3 We will simulate Scenario 3 using load steps specified in Table 5 5 AP op Xj Aw Figure 5 3 Power network system of Scenario 3 34 0 15 0 20 0 15 0 13 0 20 Table 5 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 APref 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 5 2 Design of the AGC layer for a power network using MPC The goal of the Benchmark is to design the AGC layer for the
12. instructions are described next The Controllers Cj i M are generated as follows Trucks2Dpnp is an array of N pnpmpc objects NMPC prediction horizon options sdpsettings verbose 0 Trucks2Dpnp modelTrucks2Dd createCtrlPnPMPC NMPC 2 30 options 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 29 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 13 that is free available for academic use To start a simulation we can use the function runSimTrucks2D in the following way options sdpsettings verbose 0 x0 repmat 1 0 1 0 N 1 Tsim 20 xu runSimTrucks2D Tsim x0 modelTrucks2Dc modelTrucks2Dd Trucks2Dpmp 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 1 we propose the case of N 1024 trucks and we give the results in terms of computational times 30 Chapter 5 Hycon2 Benchmark Power Network System Note for MatLab simulations in this chapter we recommend the use of CPLEX optimizer 13 5 1
13. 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 7 1 2 The createCtrIPnPMPC4lsmodel function Given an LSmodel object we can design PnPMPC controllers using 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 7 1 1 and Section 4 2 1 7 1 3 Example We create an LSmodel object of the power network system proposed in Scenario 1 in Chapter 5 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 101 10 15 kmin 8 8 8 6 1
14. proposed in 10 We illustrate the use of the class parameterizedRCI in the following example 6 3 1 examplelparameterizedRCI m matrices of the LTI system A 11 01 B 05 1 1 matrices for constraints on state and input variables cx eye 2 eye 2 dx E 22 1 52 1 cu dg dI du 3 02 3015 4 parameter of the algorithm greater than the controllability index of A B 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 control invariant set plot Other useful methods for the parameterizedRCI class are isinside Test if a point is inside the parameterized robust control invariant set R 45 x 1 1 test 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 polytope object of the set can be obtained through the attribute RP R double R RP 6 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
15. scenarios introduced in Section 5 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 PN x t 5 3a t N 1 min gt llek 2 llq lluk u llr llet N 2 ls 5 3b u t t N 1 i k t x k 1 Ax k Bu k LAPr t kKe0 N 1 5 3c x k e X ke0 N 1 5 3d u k e U kKE0O N 1 5 36 a N X 5 38 and then apply AP ef u 0 We note that the cost function depend upon x and u that are defined as xi 0 0 APr APr and ur AP The constraints X and U in 5 3d and 5 3e are obtained from constraints listed in Table 5 2 In the cost function 5 3b we set N 15 Q diag Q1 Qm and R diag R 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 35 given by APrie k Pij A0 k A6 k 5 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 Sis full MPC full we compute the symmetric positive definite matrix S and the static state feedback auxiliary control law Kauxt by maximizing the volum
16. 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 lss 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 3 3 1 examplelsubss m In this example we show how to add a subsystem to an lss object using a subss object example21lss A 0 1 k M h M 21 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 the 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 will not be saved sub will have on
17. the subsystem number i from an lss object Note that if subsystem 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 e if two subsystems 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 19 example4lss 4 remove exogenous signal number 2 modelCart modelCart removeExoSignal 2 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 the methods of the PnPM
18. toolbox 5 that allows one to plot the graph of a large scale system composed by interconnected systems Please note that MPT YALMIP and GraphViz4MatLab get automatically installed when installing PnPMPC with an option to avoid their installation in case you already have them One can check if all required toolboxes are correctly installed with the following instructions yalmiptest mpt_init 2 1 3 Installed solvers The following solvers get automatically installed when installing PhPMPC e Sedumi 6 e GLPK GNU Linear Programming Kit 7 e SDPT3 8 1 4 Installation of the PnPMPC toolbox 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 if you already have installed MPT and or YALMIP pnpmpc toolbox_init will give the possibility of skipping their installation However we guarantee that PnPMPC works correctly only with the versions of MPT and YALMIP supplied with PnPMPC There fore if you want to be 100 sure that everything works remove any previous copies of MPT and YALMIP from your disk before installing those supplied with PnPMPC Also remove any MPT and YALMIP directory from your Matlab path Remark 2 also for the solvers pnpmpc_toolbox_init will give the possibility of skipping their installation However if you want to be 100 sure that everything works remove any previous copies of the solvers from your disk before installing those
19. wW Il 145 1 1 1 17 1 03 w Il bounded disturbances as polytope object W polytope eye 2 eye 2 1111 4 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 polytope object F epsilon_mRPI A B K W epsilon 4 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 4 W zonotope 0 5 0 1 eye 2 4 F epsilon_mRPI A B K W epsilon plot the approximation of the mRPI F plot Other useful methods fot the epsilon _mRPI class are isinside Test if a point is inside the e mRPI set F ee se F isinside x doubleHK If W is a polytope object the function returns the H representation of the mRPI Moreover a polytope object of the mRPI is saved in F_epsilon F doubleHK F F_epsilon 6 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 x E X polytope eye 2 eye 2 222 2 epsilon_mRPI A B K W epsilon X 43 6 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 al
20. Hcen Kcen objlss doubleCenInputConstraint index selection Next we provide the code for adding constraints described in the code comments to the model Cart object example21lss 0 x_3 1 x_4 lt 3 modelCart modelCart addStateConstraint 2 0 1 3 4 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 L7 0 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 Vo 20u2 lt 3 modelCart modelCart addInputConstraint 2 2 3 1 u_1 2 u_2 lt 3 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 Ksu giving rise to the constraints Hsu du lt K w The method declaration is objlss objlss addDeltaUConstraint sysindex Hdeltau Kdeltau And the corresponding double method is Hdeltau Kdeltau objlss doubleDeltaUCon
21. Introduction An example of a real 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 1 within the HYCON2 project 14 We consider a PNS as composed by several power generation areas coupled through tie lines 15 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 15 Dif i Ai Bux LAP 5 Aij g 5 1 GEN where xy A0 Awi APmi AP is the state uj APref is the control input of each area AP is the local power load and M is the sets of neighboring areas i e areas directly connected 31 to DE through tie lines The matrices of system 5 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 5 2 0 000 0 Pij Aij 2H 0 0 0 Li mn IF 0 000 0 0 000 0
22. Linear Programming Kit Online Available http www gnu org software glpk glpk html K C Toh M J Todd and R H Tutuncu SDPT3 A Matlab software package for semide nite programming version 2 1 Optimization Methods and Software vol 11 pp 545 581 1999 Online Available http www math nus edu sg vmattohkc sdpt3 html S V Rakovi and D Q Mayne A simple tube controller 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 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 61 12 13 14 15 16 17 18 19 20 J B Rawlings and D Q Mayne Model Predictive Control Theory and Design Madison WI USA Nob Hill Pub 2009 IBM IBM ILOG CPLEX Online Available http www 03 ibm com ibm university academic pub page academic initiative Highly complex and networked control systems HYCON2 Network of excellence Online
23. PC 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 3 3 The subss class This class allows one to model individual subsystems and it is useful for 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 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 Many methods of this class have the same name of methods of the Iss class and work in the same way In addition many properties are
24. PnPMPC Toolbox v 0 9 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 March 2013 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 Model Predictive Control MPC schemes described in 1 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 pnpmpe_toolbox unipv it 1 1 Notations e R is the set of real number e N is the set of integers a b is the 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 posit
25. S and the ellipsoidal terminal constraint executing the function objpnpmpc objpnpmpc XFqpmax Q R options 4 1 3 The zeroTerminal method Using zeroTerminal we compute a zero terminal constraint as proposed in Chapter 2 in 12 In particular we consider Li amp 1 k vq k E k i Qt ET 1 oy un E Rlou lk vn 6 Ve q Ni 0 Xx 4177 Mi where Q Q gt 0 R and R R gt 0 e R ais 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 4 1 4 The uRH method In order to compute the constrol action uj we need to solve the optimization problem 4 3 The control action is computed executing the function u diagnostics objpnpmpc uRH x0 din xcrefin vrefin options e x0 is the initial state zj 0 a vector of dimension n x 1 for a decentralized implementation of PnPMPC controllers For a distributed implementation as in Section 5 of 1 we consider x0 as a cell array xp 0 j 0 where zga 0 is the initial state for the current subsystem and xyz 0 is the state of predecessor subsystems e xcrefin means the reference state of nominal system ref over the control horizon it can be expressed as a vector n x 1 if reference is constant over the control horizon or a matrix n x N the co
26. System 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 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 signal and a centralized input to the lss object modelCart through examplellss m With reference to 2 7 and 2 8 we want to set 0 0 0 0 0 M pe N FE Been 3 4 16 0 0 0 0 2 3 0 0 5 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 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
27. 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 reX 6 gt 42 X and e Bs 0 a x where for 6 gt 0 Bs v x e R x 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 E U C R ifVa t e 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 and w t e W C R if Va t X there exist u t U such that x t 1 X Vw t e W 6 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 16 We propose two examples 6 1 1 examplelepsilon_mRPI m The code below shows how to compute an e mRPI set for the LTI system x Ax Bu w where x R u Kx 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 42
28. art 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 number 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 4 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 maaxh The system is characterized by the matrices u 1 u 2 u N 1 u N k 2 TI dim re x 1 1 x 2 1 x N 1 1 Figure 3 2 Large scale system to model in example3
29. be performed with the MatLab command load scenario2DssZeroData PnPMPC case In the subfolders of each scenario there are MatLab files for PnPMPC simulations as scenario number scenario Full Zero De Di Data Assume we want to sim ulate Scenario 2 in a distributed fashion with zero terminal constraint We need to load MatLab file scenario2ZeroDiData This operation can be performed with the MatLab com mand load scenario2ZeroDiData e Step 3 Start a simulation from Simulink will produce the results of the control experiments sim simulatorPNS_AGC_2 5 5 Summary of results In this section we summarize the simulation results using centralized and plug and play MPC In Figure 5 4 we compare the performance of the PnP scheme with the performance of the centralized MPC controller for Scenario 1 In the control experiment step power loads APL are specified in Table 5 3 and they account for the step like changes of the control variables in Figure 5 4 We highlight that the performance of decentralized and centralized MPC are totally comparable in terms of frequency deviation Figure 5 4 a control variables Figure 5 4 b and power transfers AP e Figure 5 5 The values of performance parameter n and using different controllers are reported in Table 5 8 and Table 5 9 respectively In terms of parameter 1 plug and play controllers with decentralized and distributed online implementation are equivalent to centralized controlle
30. bles for future computations The third and fourth input argument are empty that means the user does not know convex functions gf and hf such that f gf hf Therefore the function outerApproximation computes gf and hf We highlight that the IntLab toolbox http www ti3 tuhh de rump intlab is needed Zu ZuVec Xs J H outerApproximation X f options plots figh figure 1 subplot 1 2 2 hold on plot Zu y plot ZuVec r plot xf Cisi x 241 Ba box on subplot 1 2 1 plot Xs g box on 8 3 2 example2outerApproximation m 4 definition of function f set X and sampling of f X npoints 500 56 Figure 8 5 Results of example examplelouterApproximation m X zonotope zeros 2 1 2 0 0 3 xf zeros dimension X npoints x w CL 14 0 1 x 1 0 5 x 2 exp 0 1 x 1 72 3 0 1 0 9 x 1 0 1 x 2 0 1 cos x 2 0 05 x 2 72 12 for i 1 size xf 2 Fh Il xx randpoint X xf i f xx end definition of g and h such that f g h where g and h are convex functions g x w L 0 5 x 2 0 1 0 9 x 1 0 1 cos x 2 0 05 x 2 72 J h x w L 1 0 1 x 1 exp 0 1 x 1 72 0 1 x 2 J 4 options for outerApproximation function options split ngenerators 2 options split alpha 0 5 options split max_zono 5 4 compute zonotope approximation of f X 4 Zu approximation using 1 zonotope 4 ZuVec approximation using options split max_zono zonot
31. d on APref for systems in all Scenarios are given in Table 5 2 We highlight that all parameter values are within the range of those used in Chapter 12 of 15 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 uy APrL and aj j E Nj as exogenous inputs acronym Dss In particular we note that Dss preserves the input decoupled structure of Do while D does not 5 1 1 Scenario 1 We consider four areas interconnected as in Figure 5 1 We will simulate Scenario 1 using the load steps specified in Table 5 3 32 EE io fs s w oo 04 10 0 9 0 4 0 1 os Res 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 Pa 2 P 3 Pa5 3 o3 os 1 1 or fo Table 5 2 Model parameters and constraints for systems i 1 5 AP 2 APr 3 APref 2 Awe AEref3 Po Pr P34 APra APLaA Figure 5 1 Power network system of Scenario 1 Table 5 3 Load of power APr p u for simulation in Scenario 1 APr means a step of required power hence a decrease of the
32. e of the ellipsoid de scribed by S inside the state constraints while fulfilling the matrix inequality A BKaux S A BKaux S lt Q RKauz PO he e S is block diagonal MPCdiag we compute the decentralized symmetric positive definite matrix S and the decentralized static state feedback auxiliary control law KaurT Kaur diag K Km by maximizing the volume of the ellipsoid described by S inside the state sage while fulfilling the matrix inequality A BKauz S A BKauz S lt Q RKaux K g e Zero terminal constraint M PCzero we set S 0 and X x 5 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 lt A where Tsim is the time of the simulation From 5 5 n is a weighted average of the er Q llug E ufa k IIR 5 5 BR ror between the real state and the equilibrium state and between the real input and the equilibrium input e index Tsim 1 M gt DD APriey k Ts 5 6 Toim i 1 JEN where Tsim is the time of the simulation and AFrie is the power transfer between areas i and j defined in 5 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 5 3 Control Experiments We applied the centralized MPC schemes introduced i
33. ent centralized MPC schemes for the AGC layer 5 4 Supporting MatLab files In terms of PnPMPC controllers running file makePnpmpcControllersPNS we can create controllers with different terminal regions Then we can run each simulation executing file sce nario number scenario Full Zero De Di For example if we want run the simulation for Sce nario 2 using PnPMPC controllers in a distributed fashion with zero terminal constraints we should run scenario2ZeroDi options then a file scenario2ZeroDiData is created with all results of the simulation We note that us ing PnPMPC controllers we refer to full when using a quadratic terminal region for controller Cy i M Moreover we consider discretization system by system only options is a sdpsettings object for yalmip We highlight that different performances can be achived using different solvers We also high light 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 13 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 Mat Labpoo1l before the simulations For our simulations results we refer the interested reader to 1 One can also execute all files for modeling and designed PnPMPC controllers running makeA11PNSfiles and can delete all files for PnPMPC co
34. gorithm proposed in 10 We illustrate the use of the class localControlLyapunov in the following example 6 2 1 examplellocalControlLyapunov m matrices of the LTI system A 11 01 B og d ds matrices for constraints on state and input variables cx eye 2 eye 2 dx 221 5 2 cu Ea 3 ete ds 3 3 parameter of the algorithm gt controllability index 3 x0 is a parameter of the algorithm As x0 we consider the vertices of polytope x cx x lt dx 2 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 invariant 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 1 L isinside x ulnv Given x X compute the control action u x such that x t 1 AX x zeros 2 1 44 u L uInv x double The function returns the H representation of a A contractive control invariant set More over a polytope object of the set can be obtained through the attribute LP L double L LP 6 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
35. i Vary Xj ZEN Following 9 we set in 4 1 yg t vg Olt x t a Olt where vp Olt and amp j 0 t are optimal values of the variables v j 0 and u 0 respectively 4 2 1Robust Control Invariant RCI set Consider the discrete time Linear Time Invariant LTI system x t 1 Ax t Bu 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 z t 1 X Vw t W 2Model of subsystem i without coupling terms 23 appearing in the MPC problem Ni 1 Prat min I ll n valk Vj a N 4 3a vp 0 N 1 k 0 amp 4 0 a x t i 0 Zi 4 3b ty k 1 Aink Biva k ke0 N 1 4 3c r k Xi ke0 N 1 4 3d v k E Vi ke0 N 1 4 3e Ni E 4 3f In 4 3 N N is the control horizon 4 R R4 is the stage cost Vp R R is the final cost and Xp is the terminal set Moreover from 5 3c and 5 3e the tightened constraints X and V are defined respectively as X X 0Z Vi U 08 Z 4 4 Therefore in order to design a PnPMPC controller for subsystem i we need to solve the following problem Problem P Compute nonempty RCI Z for the nominal subsystem i treating the coupling term as a disturbance Compute sets X and V We high
36. i e we split the k th column in the middle With the operator Split Z a it is possible to split the k th column of G in a desired position a i e Splitk Z a generates two sub zonotopes Z p gk 1 a g1 gka gmld Idlling lt 1 Z p gka g1 9r 1 9m d dl iny lt 1 where gx is a column of G matrix and the parameter a 0 1 Figure 8 2 shows an example of the operator Bisectx Z applied to the zonotope Z in Figure 8 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 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 51 Zn bisect Z Zn split Z gener alpha The following code shows an example where zonotope Z is splitted 2 zonotope p G bisect Z split Z k alpha figure 1 plot Z figure 2 plot Znb struct shade 0 6 figure 3 plot Zns struct shade 0 6 Figure 8 2 An example of split Bisect Z 8 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
37. ill start from an empty object and add subsystems with addSystem method The addSystem method The addSystem method adds to an Iss 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 Dias defined in 2 5 and 2 6 The arguments name and Ts 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 System Control Toolbox or a tf object transfer function object of System Control 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 system i is coupled to system j one or more blocks among Aij Bij Cy 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 j 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 add
38. in 17 We illustrate the use of the class zonotopeRCI in the following example 6 4 1 examplelzonotopeRCI m matrices of the LTI system r T14 1 d C 0 11 matrices for constraints on state and input variables XS 4 ax ex x lt dx U u cu u lt du Os 10 9 0 42 385 3 3 3 202 13 11 3 parameter of the algorithm greater than the number of 5 matrices of the zonotope for the description of set W w f Ed d inf lt 10 01 0 3 0 4 create the zonotopic robust control invariant set 46 qa 1 qb i qp i Z zonotopeRCI A B cx dx cu du f E k qa qb qp plot the zonotopic 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 ulInv x 47 Chapter 7 PnPMPC and others toolboxes 7 1 Integration with the toolbox WIDE 18 7 1 1 The Ismodel2lss 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
39. ive definite e A diag A411 Ass is the block diagonal matrix An 0 0 0 0 1 1 0 0 Ags Zi Ul Yi dij U cen are column vectors of suitable dimensions e x Ti s means that x is composed by column vectors aj i 1 s stacked in a single column i e zi 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 Ve 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 The PnPMPC toolbox requires the following toolboxes to be installed Control System Toolbox 2 which implements the class ss state space We also use the function c2d for time discretization MPT toolbox 3 which implements the polytope class YALMIP toolbox 4 which include the optimization function solvesdp that is called for solving optimization problems e GraphViz4MatLab
40. 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 plot Zu y plot ZuVec r plot zf 1 xf 2 1 bet box on subplot 1 2 1 plot Xs g box on 59 N Figure 8 7 Results of example example3outerApproximation m 60 Bibliography 1 S Riverso M Farina and G Ferrari Trecate Plug and Play Model Predictive Control 10 11 based on robust control invariant sets Dipartimento di Ingegneria Industriale e dell Informazione Universit degli Studi di Pavia Pavia Italy Tech Rep 2012 Online Available arXiv 1210 6927 Mathworks Control System Toolbox for Matlab M Kvasnica P Grieder and M Baoti Multi Parametric Toolbox MPT 2004 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 Computer 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 http code google com p graphviz4matlab J F Sturm Using SeDuMi 1 02 A Matlab toolbox for optimization over symmetric cones Optimization Methods and Software vol 11 no 1 4 pp 625 653 1999 GNU GLPK GNU
41. light that Problem P can be solved using efficient procedure proposed in 10 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 we declare that it is impossible to design a PnPMPC controller for subsystem 7 The interested reader is refer to 1 4 1 1 The pnpmpc method Assume subss stores the model of subsystem i given by 2 3 and 2 4 Controller Cj is created through objpnpmpc pnpmpc subss N Xj k options where N is the receding horizon Xj is a cell array Hx Kxj such that neighbors states 21 j N are in the set Xj xj Hxjxj lt Kz kis a parameter of Algorithm 1 in 1 and must be gt of the controllability index of the pair A Bi One can also set a cell array kmin kmax in order to try different k options is used to set the sdpsettings for yalmip see yalmip manual for details 4 We propose two methods to compute the terminal region and the terminal penalty 4 1 2 The XFqpmax method Using XFqpmax we compute a quadratic terminal region using procedures proposed in 11 In particular we consider 24 Vi a Ni By Ni 2657 N Sla Na ST N 8 ER 8 58 lt 1 where Q Q gt 0 e R R R gt 0 RX and S S gt 0 e Ri ary k and visi k are the state and input reference trajectories for tracking capabilities We can design the terminal penalty
42. lss 0 1 0 1 An E k 1 na ANN k N 1 ud MI MO M N MN Vi 2 N 1 0 0 Ai kara _aG y n gt Ba 100 3 5 M i M i M i 0 0 0 0 Atena k i 1 _ h i 1 Agfa E k i 1 aal M i 1 M i 1 M i M i Ca 1 o Dadi The MatLab code for the model creation is when N 1000 is minM 1 maxM 10 mink 0 1 maxk 0 9 minh 0 1 maxh 10 M minM rand 1 N maxM minM 15 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 input constraints to an lss model The addition of con straints on output exogenous signal and centralized input can be done in a similar way Relevant methods a
43. lumns 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 vrefinin means the reference input of nominal system Vref over the control horizon it can be expressed as a vector m x 1 if reference is 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 25 As output we have e u control input value e diagnostics information about the optimization problems 4 2 Design of a PnPMPC controller for a large scale system PnPMPC controllers Cj for each subsystem in an Iss object are created with the method createCtr1PnPMPC described next 4 2 1 The createCtrIPnPMPC method This method acts on an Iss object and it is called as follows ctrl createCtr1lPnPMPC objlss N k option where N options have the same meaning as in the pnpmpc method The input k is a cell array of two elements kmin kmax both kmin and kmax can be scalars and then they will be used for all execution of pnpmpc function for the i subsystem or they can be vectors of dimension lxnumSys where kmin i and kmax i are used for the execution of pnpmpc function for the i th s
44. ly 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 22 Chapter 4 The pnpmpc class 4 1 Introduction In this chapter we introduce the pnpmpc class in order to design a PnPMPC controller for a single subsystem PnPMPC controllers are described in 1 In 1 the authors propose a PnP design procedure hinging on the notion of tube MPC 9 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 i M is given by 2 3 and 2 4 and we consider that each subsystem 7 is input decoupled i e Bi 0 i j Our control design procedure will hinge on the following assumptions Assumption 1 The matrix pairs Au Bi Vi e M are controllable Assumption 2 Constraints X and U i M are compact and convex polytopes containing the origin in their nonempty interior Let Z be an RCTI set for the nominal subsystem i where the coupling term jcy AijX is treated as a disturbance We design the following controller for subsystem i Ci Up va Rila Fy 4 1 where K Zi U is any feedback control law guaranteeing ty Zi gt xt t Z
45. me and method is the discretization technique All the discretization techniques of the function c2d of the Control System Toolbox which converts continuous time sys tems to discrete time models are implemented like zoh zero order hold or Tustin seehelp c2d for more help The c2d discretization methods are not always a good choice for LSS because for exact discretization 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 That means that each subsystem in 2 5 and 2 6 is discretized considering states x1 j 7 i as exogenous inputs We can save computational time because we discretize sequentially subsystems that usually have low dimension and we do not create new couplings among subsystems For using this technique set subsystem as method Other utilities There are some other object s properties that we have not explained yet e coupling is a matrix with size numSys x numSys composed by Boolean elements element i j is 1 only if subsystem i is dynamically coupled or input coupled to subsystem j i e i j ee
46. n the previous section to 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 5 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 36 it We highlight that violation of this assumption can impact considerably on the index n In all experiments we use Tsim 100 In Table 5 8 and 5 9 the values of the performance parameters 7 and respectively are reported for each control experiment See Scenario 2 Scenario 3 gt d gt Pa D De cali EPC Fall 00245 00218 00348 00347 00510 00511 PM PCaiag 0 0219 0 0249 00846 0 0547 00510 0 051 M PCzero 00249 00248 0 08486 00347 00510 0 051 Table 5 6 Values of the performance parameter n using different centralized MPC schemes for the AGC layer __ Scenario 1 Scenario 2 Scenario 3 PE D Ps D Dss D Dss PC Tal 00030 0 0020 0 0063 0 0060 10060 0 005 M PCdiag 00030 0 0030 0 0063 00061 0 060 0 0055 MPC zero 0 0030 0 0025 0 0068 0 0059 00059 0 0055 Table 5 7 Values of the performance parameter using differ
47. n zonotope and polytope sets Minkowski sum Pontryagin difference intersection N union U and relational operators C C D 2 and Zonotope arrays are managed as polytope arrays in MPT2 In the following sections we introduce some useful operators for zonotope sets Most of the definition are from 19 8 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 do gi ay 124 3141 zonotope p G plot Z The zonotope is shown in Figure 8 1 In oder to access the G representation of zonotope Z one can use the command doublePG p G doublePG Z In oder to access the H representation of zonotope Z one can use the command doubleHK H K Zo doubleHK Z 50 Figure 8 1 Zonotope Z 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 8 2 Functions for zonotope objects 8 2 1 Bisection and Split operators The operator Bisectx generates two sub zonotopes from one zonotope In particular given a zonotope Z the operator Bisect Z generates the two sub zonotopes Jk Ik Z p 3 4 gt 91 Imld dl ling lt 1 Jk Jk where gx is the k th column of G and is the biggest generator
48. nal Ma C N is the set of exogenous signals that act on system i and Mij e R 1 Nij e RPi We further enhance our model by introducing centralized inputs ucen j 80 that the dynamics of subsystem becomes ct At Buug DO Agr Bau DI Myd DI Boensijtcent 2 5 JEN JEN ai JENu i yp Cup Daug DO Cat Dyu DI Nagdy DI Deen ijttcen j 2 6 SEN JEN di IEN uri where Nu C N is the set of centralized inputs Ucen j gt J R that act on system and Been ij R and Deen RP The difference between uj and ucenfi is that uj is a local input i e the output of a local regulator for subsystem 7 while ucen j is a global input that cannot be computed locally From models 2 5 and 2 6 the collective dynamics of the resulting LSS is xt A Bu Md BeenUcen 2 7 y Cx Du Nd DeenUcen 2 8 where x zu Tj sees X s ER 2 9 is the overall state and n J jem Ni u uri Uy ey U s Ee R 2 10 is the overall input and m J jem Mi y Un Y gt Yis R 2 11 is the overall output and p Doc yy Pi deR 2 12 collects the exogenous signals acting on the overall system and Ucen R 2 13 collects the centralized inputs acting on the overall system Moreover one has Au una Ais B kae Bis Agi sea Agg By sae Bas All other matrices in 2 7 and 2 8 have a similar block structure Moreover has C R D RP Me R 4 N
49. ntrollers running clearAl1PNSfiles 37 For each control experiment we provide a file mat of the simulation that contains e Iss object of the continuos linear system pnsCn where n is the number of the scenario e parameters of the control experiment Tsim and deltaPL where deltaPL corresponds to APz e the results of the control experiment x deltaPref n and where deltaPref 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 5 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
50. o mass i in the horizontal respectively vertical 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 2j j o0 lt 1 5 j 1 3 zugllo lt 0 8 M 1 2 4 and with the input constraints upjllo 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 x j 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 28 ptions 40 30 20 10 0 10 20 30 40 50 Figure 4 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 makeModelTrucks2D as follows N 4 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 lss objectes We can design local controllers executing the function options sdpsettings verbose 0 Trucks2Dpnp makePnpmpcTrucks2D modelTrucks2Dd options whose
51. olbox 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 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 LTI 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 up R is the local input at time k Moreover we have A R Ai R i By E 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
52. opes 4 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 57 hold on plot Zu y plot ZuVec r plot xf 1 xf 2 b box on subplot 1 2 1 plot Xs g box on Figure 8 6 Results of example example2outerApproximation m 8 3 3 example3outerApproximation m definition of function f set X and sampling of f X differently from example2outerApproximation we want to split only along x 1 therefore x 2 is substituted by w 1 npoints 500 X1 zonotope 0 2 W 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 1 for i 1 size xf 2 xx randpoint X xfi i f xx 1 xx 2 end 58 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 1 0 1 x 1 exp 0 1 x 1 72 O 1 w 1 195 4 options for outerApproximation function options split ngenerators 2 options split alpha 0 5 options split max_zono 5 4 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
53. ough a nonlinear function using DC programming which is based on DC functions 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 a and h x are convex functions In order to compute a tighter outer approximation in 19 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 propose three examples that can be found in the PnPMPC toolbox 1The function outerApproximat ion needs Symbolic Math Toolbox and Partial Differential Equation Toolbox 55 8 3 1 examplelouterApproximation m 4 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 1 for i 1 size xf 2 xx randpoint X xf i f xx end 4 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 i zonotope ZuVec approximation using options split max_zono zonotopes Xs zonotope X splitted J jacobian and H hessian are outputs varia
54. plemented in PnPMPC toolbox in the following instruction Zo reduce Z ngen The following code shows an example where the number of generators of zonotope Z is reduced C1 1 1242 45 31456 zonotope p G ngen 3 Zo reduce Z ngen figure 1 hold on plot Zo r plot Zz b 20 10 5b 10F Figure 8 4 The reduction operator redngen Z is applied to the green zonotope yielding a more conservative approximation red zonotope 8 2 4 The support function of polytope zonotope sets The support function of a polyhedral set is defined as 8 1 sup max cx 8 1 For polytope and zonotope objects we developed the function supportFunction One can compute sup using the following instruction sup supportFunction X c where X is a polytope or zonotope object 8 2 5 Zonotope Polytope e Polytope P is a zonotope bool iszonotope P e Polytope P to Zonotope Z Z polytope2zonotope P e Zonotope Z to Polytope P P zonotope2polytope Z 8 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 20 19 In 20 an algorithm to compute zonotope outer approximations for nonlinear systems was proposed The authors suggest to create an image of a zonotope thr
55. pmpcTrucks1D or modelTrucks1Dd createCtr1PnPMPC 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 13 that is free available for academic use To start a simulation we can use the function runSimTrucks1D in the following way options sdpsettings verbose 0 solver sdpt3 x0 repmat 1 0 N 1 27 Tsim 20 x u runSimTrucks1iD Tsim x0 modelTrucksiDc modelTrucksiDd TrucksiDpnp This function computes the control action using the method uRH of the pnpmpc class We used sdpt3 solver that is installed with the toolbox but we really recommend the use of CPLEX 4 3 2 Example 2D We consider a large scale system composed by N masses coupled as in Figure 4 3 2 for the case of N 1024 where the four edges connected to a point correspond to springs and dampers arranged as in Figure 4 2 Hereafter we assume that N 2 for some z N z gt 0 Each Figure 4 2 Array of masses details of interconnections mass i M 1 N is a subsystem with state variables x 2 i 1 2 2 i 3 Tji 4 and input up Uti 1 Ufs 2 where 1 and 2 3 are the displacements of mass i with respect to a given equilibrium position on the plane equilibria lie on a regular grid x 2 and ya are the horizontal and vertical velocity of the mass i respectively and 100u 1 respectively 100u 9 is the force applied t
56. r however the performance of PnP DeMPC are such that each area can absorb the local loads by producing more power locally AP 6 i instead of receiving power from predecessor areas for this reason PnP DiMPC has performance more similar to centralized MPC The values of performance parameter n and using different controllers for Scenario 2 and 3 are reported in Table 5 8 and Table 5 9 respectively We refer the interested reader to 1 for an in depth analysis of the results for all scenarios seen Scenario Scenario Paio nt Cem MPOding MPCsero MPCiiog MPCzer amie 0020 009 oo omar ooo 001 De PuPMPG 2 81 281 1502 490 _ 7 65 765 DiPaPMPC 3 61 1301 231 231 2 15 215 Table 5 8 Value of the performance parameter n for centralized MPC first line and percentage change using decentralized and distributed MPC schemes for the AGC layer Best values for PnP controllers are in bold sceso Scanio Senarios __ C PCa uP Caer MPCdiny MPC MP Ching MPC Cene 000 00029 0006 0001 osso 0008 De PaPMPC 26 667 24 14 31 25 27 08 42 85 58 097 DiPNPMPC 000 1344 862 517 Tuam 357 Table 5 9 Value of the performance parameter for centralized MPC first line and percentage change using decentralized and distributed MPC schemes for the AGC layer Best values for PnP controllers
57. randomly chosen in the interval 1 1 5 We obtain models Xj 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 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 In order to create the large scale system composed by N masses we can use the function makeModelTrucks1D as follows N 4 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 options sdpsettings verbose 0 TrucksiDpnp makePnpmpcTrucks1D modelTrucksiDd options whose instructions are described next The Controllers Cj i M are generated a follows Trucks1Dpnp is an array of N pnpmpc objects NMPC 20 prediction horizon options sdpsettings verbose 0 TrucksiDpnp modelTrucksiDd createCtrlPnPMPC NMPC 1 30 options Then we add a zero terminal constraint to each controller using the following instructions Q 10 eye 2 R 1 parfor i 1 modelTrucksiDd numSys TrucksiDpnp i Trucks1iDpnp i zeroTerminal Q R end We note that the functions makePn
58. re 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 matrices 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 inputs outputs exogenous signals and centralized 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 MPT toolbox 3 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 inde
59. 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 predecessors to subsystem i is N j i j x 1 SubSystem 1 SubSystem 3 x 3 x 2 SubSystem 2 SubSystem 4 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 successors of subsystem 7 i e the set of subsystems influenced by it Definition 6 The set of successors to system i is Si j j i E E For example from Figure 2 1 one has S2 4 since A42 0 The output of system is given by Yi Curly Diuji 5 Citi Diju 2 2 SEN where yp RP Cu RP Dj RE Ci RP and Di e RPM Definition 7 Two systems 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 Diu DI Nidy 2 4 JEN JENa where dy R is an exogenous sig
60. straint index selection 3 2 5 example5lss 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 system 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 B C D respectively If how is not passed blocks i j in all matrices A B C D will be removed See help removeCoupling for more information The removeExoSignal method objlss objlss removeExoSignal i This method allows one to remove the exogenous signal dj It deletes the column i of matrices M and N and also the constraints related to this exogenous signals if present Note that a coupling 18 constraint related to exogenous signals and j after the removal of the signal becomes a local constraint of signal j For example let dj 2d a 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
61. supplied with PnPMPC Also remove any solver directory from your MatLab path Remark 3 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 5 Directories The PnPMPC toolbox consists of the following directories e pnpmpc_toolbox epsilon_mRPI contains methods for computing outer approximation of minimal robust positively invariant sets e pnpmpc_toolbox localControlLyapunov contains methods for computing control invariant sets 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 pnpmpc contains methods for design of PnPMPC controllers e pnpmpc_toolbox polytope contains functions of MPT toolbox improved to use zono tope sets 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 docs contains html documentation pnpmpc_toolbox examples contains examples demonstrating the functionalities of PaPMPC toolbox pnpmpc_toolbox extra contains extra useful functions pnpmpc _toolbox solvers contains solvers provided with PnPMPC toolbox pnpmpc_toolbox to
62. tputs of every subsystem as in 2 9 2 10 and 2 11 numExo numICen are respectively nq in 2 12 and n 2 13 all constraint matrices H and K appearing in 2 21 Hdeltau Kdeltau that will be explained in Section 3 2 4 coupling name that will be explained in Section 3 1 namex nameu namey named nameucen in order to assign names to variables see Section 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 of MatLab sparse matrices by default This usually allows one to save a considerable amount of memory To visualize the matrices in the full format is enough to use the MatLab command full 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 is no exogenous signal or centralized input 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 ti
63. ubsystem This method extracts a subsystem calls the method pnpmpc and creates the controller The output ctrl is an array of dimension 1xnumSys and in the position 7 it is stored the controller for subsystem i 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 4 3 Examples As an example of real application in Chapter 5 we describe PnPMPC controllers for power network systems and run simulations using methods proposed in this chapter In this section we propose simpler example 4 3 1 Example 1D We consider a large scale system composed by N masses coupled as in Figure 4 1 Each mass u 1 u 2 u N 1 u N k 2 k N 2 Imm am x 1 1 x 2 1 x N 1 1 x N 1 Figure 4 1 Array of masses in 1D i M 1 N is a subsystem with state variables xj u 2 and input up Uj 26 where 1 is the displacements of mass 7 with respect to a given equilibrium position x 2 is the horizontal velocity of the mass and 1004 1 is the force applied to mass 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 x 1 lo lt 1 5 p a Iloo lt 0 8 i M and with the input constraints lu lt Bi where 8 have been
64. ure 3 1 The model is 1 1 x 2 1 Figure 3 1 System modeled in examplellss m Ti T 2 wii UD B 3 1 where zp is composed by 1 1 the displacement of first mass with respect to a given equilibrium position and 1 9 the velocity of the first mass Similarly x is composed by 22 1 and 22 3 Moreover u and uf are forces applied to the first and second mass respectively We also have 0 1 h k M M Au A19 Agi A22 A Ay A22 0 0 Aw A2 ie Hl 3 2 M M B diag B11 B22 By Boo do 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 the sizes of blocks Aij and Bij in matrices A and B Cand D are the matrices in 2 8 pi is a vector storing the sizes of blocks Cij and Dij in matrices C and D 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 w
65. x 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 returned 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 polytope i e an object of polytope class see the documentation of the MPT toolbox for help of system specified in index giving rise to a state constraint Also if index contains more than one index a polytope 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

PnPMPC Toolbox v. 0.9 - User manual

Contents

Download Pdf Manuals

Related Search

Related Contents