A Homogeneous And Self-Dual Interior-Point Linear Programming Algorithm .
SUBMITTED FOR IEEE TRANSACTIONS ON AUTOMATIC CONTROL, AUGUST 20, 20141A Homogeneous and Self-Dual Interior-PointLinear Programming Algorithm for EconomicModel Predictive ControlLeo Emil Sokoler, Gianluca Frison, Anders Skajaa, Rasmus Halvgaard, John Bagterp JørgensenAbstract—We develop an efficient homogeneous and self-dualinterior-point method (IPM) for the linear programs arising ineconomic model predictive control (MPC) of constrained linearsystems with linear objective functions. The algorithm is basedon a Riccati iteration procedure, which is adapted to the linearsystem of equations solved in homogeneous and self-dual IPMs.Fast convergence is further achieved by means of a recentwarm-start strategy for homogeneous and self-dual IPMs. Weimplement the algorithm in MATLAB and C. Its performanceis tested using a conceptual power management case study.Closed loop simulations show that 1) the proposed algorithmis significantly faster than several state-of-the-art IPMs basedon sparse linear algebra, and 2) warm-start reduces the averagenumber of iterations by 35-40%.Index Terms—Optimization algorithms, Linear programmingalgorithms, Predictive control for linear systems, Riccati iterations, Energy systemsstructured LP. Section III introduces the homogeneous andself-dual model. This section also provides a Riccati iterationprocedure for solving the structured linear system that occursin every iteration of an IPM for solving the LP. SectionIV discusses warm-start. Section V compares a MATLABand C implementation of the proposed algorithm denotedLPempc to several state-of-the-art IPMs using a simple powermanagement case study1 . Concluding remarks are given inSection VI.II. O PTIMAL C ONTROL P ROBLEMWe consider linear state space systems in the formxk 1 Axk Buk Edk ,yk Cy xk ek ,The authors are with the Department of Applied Mathematics and ComputerScience, Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark.L. E. Sokoler and A. Skajaa are also affiliated with DONG Energy, DK-2820Gentofte, Denmark (e-mail: {leso, giaf, andsk, rhal, jbjo}@dtu.dk)(1a)ek N (0, Re ),(1b)zk Cz xk ,I. I NTRODUCTIONIn economic MPC of linear systems with a linear objectivefunction and linear constraints, the constrained optimal controlproblem can be posed as a linear program (LP). As theoptimization problem is solved online, the performance andreliability of the optimization algorithm solving the LP isimportant. In this paper, we develop a homogeneous and selfdual variant of Mehrotra’s predictor-corrector method [1], [2]for economic MPC that combines the following performanceimproving components: Riccati Iteration Procedure: We speed-up the most timeconsuming numerical operations using a Riccati iterationprocedure. Warm-Start: We implement a warm-start strategy for thehomogeneous and self-dual IPMs. [3] reports that thismethod reduces the number of iterations by 30-75% basedon the NETLIB collection of test problems.Riccati-based IPMs have been developed for set-point basedMPC with an 2 -penalty [4]–[6] and with an 1 -penalty [7]. ARiccati iteration procedure has not previously been combinedwith the homogeneous and self-dual model. The homogeneousand self-dual model has become widely adopted by state-ofthe-art IPMs for linear and conic programming. This paperis organized as follows. Section II formulates the control lawassociated with economic MPC as the solution to a highlydk N (0, Rd ),(1c)where x0 N (x̂0 , P0 ). (A, B, Cy , Cz , E) are the state spacematrices, xk Rnx is the state vector, uk Rnu is the inputvector, yk Rny is the measurement vector, zk Rnz isthe output vector, dk is the process noise vector and ek isthe measurement noise vector. We use bold letters to denotestochastic variables.In this paper, the Economic MPC optimal control problemis defined to optimize the control actions of (1) with respect toa linear economic objective function, input limits, input-ratelimits and soft output limits. Evaluation of this control lawrequires the solution to the LPXTmin.pTk j uk j qk j 1ρk j 1 ,(2a)u,x̂,ẑ,ρj N0s.t.x̂k j 1 k Ax̂k j k Buk j ,j N0 ,ẑk j k Cz x̂k j k ,j N1 ,(2c)uk j uk j uk j ,j N0 ,(2d) uk j uk j uk j 1 uk j ,j N0 ,(2e)z k j ρk j ẑk j k z k j ρk j ,j N1 ,(2f)ρk j 0,j N1 ,(2g)(2b)where Ni : {0 i, 1 i, . . . , N 1 i}, with N beingthe length of the prediction and control horizon. The problemdata are the state-space matrices, (A, B, Cz ), the filtered1 Software and implementation details are available via http://www2.imm.dtu.dk/ jbjo/publications
SUBMITTED FOR IEEE TRANSACTIONS ON AUTOMATIC CONTROL, AUGUST 20, 2014estimate, x̂k k , the input limits, (uk j , uk j ), the inputrate limits, ( uk j , uk j ), the output limits, (z k j , z k j ),the input prices, pk j , and the price for violating the output constraints, qk j . Note that for compact notation, theoptimization variables in (2) are written T as u, x̂, ẑ and ρ,where u uTk uTk 1 · · · uTk N 1 , and similarly forx̂, ẑ, and ρ. The filtered estimate, x̂k k : E[xk Yk ], isthe conditionalexpectation of T xk given the observationsYk : y0T y1T y2T . . . ykT . We obtain this value usingthe Kalman filter. By augmenting the state-space system suchthat A 0x̂kBA : , x̂k : , B : ,0 0uk 1I EE : , Cz : Cz 0 , Cy : Cy 0 ,0we can express (2e) asj N0 , uk j uk j Dxk uk j , in which D : 0 I . This formulation simplifies latercomputations considerably. To keep the notation simple weassume that k 0 and write x̂j : x̂0 j 0 for conditionalexpressions. The problem data is collected in the structures g,F , b, H and c, and (2) is put into the formmin. {g T t F t b, Ht s c, s 0}.(3)t,sAs an example, consider the case for N 2 Tt : uT0 x̂T1 ρT1 uT1 x̂T2 ρT2 , Tg : pT0 0 q1T pT1 0 q2T ,and F b : H c : B0 IA000B0 I00 Ax̂00u0I00000000I00u1 I00000 u0000 I00 u1I00000 ũ00 D 0I00 u1 I00000 u0 0D0 I00 u10Cz I 000z10000Cz Iz20 Cz I 000 z 10000 Cz I z 200 I 000000000 I0 , , where ũ0 : u0 Dx̂0 and u0 : u0 Dx̂0 . The problem (2) can thus be posed as a highlystructured LP withn : N (nu nx nz ) variables, mE : N nx equalityconstraints, and mI : N (4nu 3nz ) inequality constraints.Note that we have eliminated ẑj from the optimization problemusing the linear relation (2c).2III. H OMOGENEOUS AND S ELF -D UAL I NTERIOR -P OINTM ETHODThe dual of the LP (3) ismax { bT v cT w F T v H T w g, w 0},v,w(4)where v RmE and w RmI are dual variables corresponding to the Lagrange multipliers for the equality constraints andthe inequality constraint of (3), respectively. We assume thatF has full row rank. This is always the case for the problem(3).Homogeneous and self-dual IPMs solve (3) and (4) indirectly. The idea is to construct a related LP with propertiesthat are useful for IPMs. Aside from an inherent ability todetect infeasibility, recent advances show that IPMs based onthe homogeneous and self-dual model can be warm-startedefficiently [3]. We refer to [8]–[10] for proofs and details.Introduce a new set of optimization variables (t̃, ṽ, w̃, s̃),and the additional scalar variables (τ̃, κ̃). Then the self-dualand homogeneous problem associated with (3) and (4) may bestated as the linear feasibility problemfind t̃, ṽ, w̃, s̃, τ̃, κ̃,TTs.t. F ṽ H w̃ gτ̃ 0,(5a)(5b)bτ̃ F t̃ 0,(5c)cτ̃ H t̃ s̃ 0,(5d) g T t̃ bT ṽ cT w̃ κ̃ 0,(5e)(w̃, s̃, τ̃, κ̃) 0,(5f)Proposition 1 shows that the solution to (3) and (4) can beobtained by solving (5).Proposition 1: The linear feasibility problem (5) alwayshas a strict complimentary solution (t̃ , ṽ , w̃ , s̃ , τ̃ , κ̃ )satisfying s̃ j w̃j 0 for j 1, 2, . . . , mI and τ̃ κ̃ 0.For such a solution, one of the following conditions hold1) τ̃ 0, κ̃ 0: The scaled solution (t , v , w , s ) (t̃ , ṽ , w̃ , s̃ )/τ̃ is a primal-dual optimal solution to(3) and (4).2) τ̃ 0, κ̃ 0: The problem (3) is infeasible or unbounded; either bT ṽ cT w̃ 0 (primal infeasible),or g T t̃ 0 (dual infeasible).Proof: See [10], [11].A. Riccati Iteration Procedure for Economic MPCLet k denote a particular iteration number. The necessaryand sufficient optimality conditions for (5) are (w̃, s̃, κ̃, τ̃ ) 0and F T ṽ H T w̃ gτ̃0 0 bτ̃ Ft̃ 0 T cτ̃ T H t̃ Ts̃ (6) g t̃ b ṽ c w̃ κ̃ 0 , 0 W̃ S̃1mI0τ̃ κ̃W̃ is a diagonal matrix with the elements of w̃ in its diagonal,and similarly for S̃. Moreover, 1mI is the column vector ofall ones of size mI . The main computational bottleneck in
SUBMITTED FOR IEEE TRANSACTIONS ON AUTOMATIC CONTROL, AUGUST 20, 2014finding a point that satisfies the optimality conditions usinga homogeneous and self-dual IPM is solving a linear systemof equations (in each iteration of the IPM) [1], [2]. For anarbitrary right hand side, we can solve this linear system ofequations asF T ṽ H T w̃ g τ̃ r1 ,(7a)b τ̃ F t̃ r2 ,(7b)c τ̃ H t̃ s̃ r3 ,(7c)g T t̃ bT ṽ cT w̃ κ̃ r4 ,(7d)W̃ s̃ S̃ w̃ r5 ,(7e)κ̃ τ̃ τ̃ κ̃ r6 .(7f)The system (7) is different from the system solved in standardIPMs. The Riccati iteration procedures proposed in [4]–[6]can therefore not be applied directly to solve the system (7).Proposition 2 shows that a solution to (7) can be computedby solving a reduced linear system of equations (8) and anumber of computationally inexpensive operations. Proposition 3 shows that the reduced system can be solved using aRiccati iteration procedure that scales as O(N (nu nx nz )3 ).The procedure eliminates variables from the reduced linearsystem (9), such that it can be expressed in a form (11) thatoften occurs in optimal control problems. Riccati iterationssolve the resulting system efficiently by exploiting its structure[4]–[6]. The complexity of solving the system (7) directlyusing sparse linear algebra routines is linear to quadratic inN , while a general purpose solver using dense linear algebrascales cubically [12]. If the number of states, nx , is largecompared to the number of inputs, nu , condensing methodsare more efficient than Riccati-based solvers [5]. Note that itis convenient to have an algorithm that scales linearly in theprediction horizon, N , as stability of MPC schemes often maybe achieved by selecting a sufficiently large value of N [13],[14].Reference [15] provides details on a predictor-corrector IPM[1], [2], LPempc, that utilizes a C implementation (with callsto BLAS and LAPACK) of the numerical procedure presentedin Proposition 2 and Proposition 3. To speed-up the numericalcomputations and reduce the storage requirements, LPempchandles operations involving the structured matrices F and Has specialized linear algebra routines.Proposition 2: The solution to (7) can be obtained bysolving HT0 FTf1 h1r1 g F 00 f2 h2 r2 b ,(8)f3 h3r3 c H 0 W̃ 1 S̃and subsequent computation ofr6 τ̃ (g T f1 bT f2 cT f3 ), τ̃ κ̃ τ̃ (g T h1 bT h2 cT h3 ) t̃ f1 h1 τ̃,3where r3 : r3 W̃ 1 r5 and r6 : r6 τ̃ r4 .Proof: See [11].Proposition 3: System (8) can be solved in O(N (nu nx nz )3 ) operations using a Riccati iteration procedure.Proof: For a single arbitrary right hand side, we maywrite the system (8) as 0 FTHTrD t̃ F 00 ṽ rE .(9)rI w̃ H 0 W̃ 1 S̃Denote the Lagrange multipliers associated with the inequalityconstraints (2d) and (2g) as η, δ, υ, ω, γ, ζ and ξ where TT,η : η0T η1T . . . ηN 1and similarly for δ, υ, γ, ζ and ξ. The multipliers (η, δ) areassociated with the input limits (2d), (υ, ω) are associatedwith the input-rate limits (2e), (γ, ζ) are associated with theoutput limits (2f), and ξ is associated with the non-negativeconstraints (2g). The system variables are written in the form T t̃ uT0 x̂T1 ρT1 . . . uTN 1 x̂TN ρTN , TT ṽ ṽ0T ṽ1T . . . ṽN, 1 T T w̃ η δ T υ T ω T γ T ζ T ξ T .Accordingly, we partition the right hand side such that T TTTTTTrx,1rw,1. . . ru,N,rD ru,0 1 rx,N rw,N T TTTrE Rv,0 Rv,1 . . . Rv,N 1 , T rI rηT rδT rυT rωT rγT rζT rξT ,and write the diagonal matrix W̃ 1 S̃ in terms of diagonalsubmatrices W̃ 1 S̃ diag ΣTη , ΣTδ , ΣTυ , ΣTω , ΣTγ , ΣTζ , ΣTξ .The linear system of equations (9) can now be stated in theform ηi δi υi ωi B T ṽi ru,i , ṽi ui Ση,i ηi rη,i ,i N0 , ui Σδ,i δi rδ,i ,i N0 , ui D x̂i Συ,i υi rυ,i ,i Ñ0 , ui D x̂i Σω,i ωi rω,i ,i Ñ0 , x̂i 1 A x̂i B ui Rv,i ,i Ñ0 , ρi Cz x̂i Σγ,i γi rγ,i ,i N1 , ρi Cz x̂i Σζ,i ζi rζ,i ,i N1 , ρi Σξ,i ξi rξ,i ,i N1 , γi ζi ξi rw,i ,i N1 ,CzT ( γi 1T ζi 1 ) A ṽi DT (ωi υi ) rx,i ,with Ñ0 : N0 \ {0} and the special cases ṽ f2 h2 τ̃, u0 Συ,0 υ0 rυ,0 , w̃ f3 h3 τ̃, u0 Σω,0 ω0 rω,0 , κ̃ g T t̃ bT ṽ cT w̃ r4 , s̃ W̃ 1(rC S̃ w̃),i N0 , ṽN 1 x̂1 B u0TCz ( γN ζN ) Rv,0 , rx,N .i Ñ0 ,
SUBMITTED FOR IEEE TRANSACTIONS ON AUTOMATIC CONTROL, AUGUST 20, 2014The Lagrange multipliers associated with the inequality constrains η, δ, υ, ω, γ, ζ and ξ are eliminatedfrom the equations above. This is computationally inexpensiveas the matrices to be inverted in the process are all diagonal.The reduced set of equations areB T ṽ0 U0 u0 Ru,0 ,(10a)B T ṽi Ui ui Gi x̂i Ru,i ,i Ñ0 , (10b) x̂1 B u0 Rv,0 ,(10c) x̂i 1 A x̂i B ui Rv,i ,i Ñ0 , (10d)Wi ρi MiT x̂i Rw,i ,i N1 , (10e) ṽi 1 Mi ρi X̄i x̂i GTi ui AT ṽi R̄x,i ,i Ñ0 , ṽN 1 MN ρN X̄N x̂N R̄x,N ,(10f)(10g)where we have definedUi : Gi : 1 (Σ 1ω,i Συ,i )D, 1 1Σ 1ζ,i Σξ,i Σγ,i , 1CzT (Σ 1ζ,i Σγ,i ), 1CzT (Σ 1ζ,i Συ,i )CzWi : Mi : X̄i : X̄N : CzT (Σ 1ζ,N This system can be solved by a Riccati iteration procedurein O(N (nu nx nz )3 ) operations [4], [6], [7].IV. WARM -S TARTWe apply the warm-start strategy of [3] to pick an initialpoint for the homogeneous and self-dual IPM. The methodcombines a guess of the solution (candidate point), (t̄, v̄, w̄, s̄),with the standard starting point (0, 0, 1mI , 1mI , 1, 1). Animportant feature of the homogeneous and self-dual modelis that the the standard starting point is perfectly centralizedwith respect to the central path [3]. This makes warm-startwork well for IPMs based on the homogeneous and self-dualmodel. The initial point in [3] is defined asw̃0 λw̄ (1 λ)1mI ,t̃0 λt̄,0Σ 1η,i Σ 1δ,i4Σ 1ω,i τ̃ 1,Σ 1υ,i ,i N0 ,i Ñ0 ,i N1 ,i N1 , DT(Σ 1γ,i Σ 1ω,i )D,i Ñ0 ,Σ 1υ,N )Cz ,andRu,i : ru,i r̄δ,i r̄ω,i r̄η,i r̄υ,i ,i N0 ,Rv,i : Rv,i ,i N0 ,Rw,i : rw,i 1 r̄ζ,i 1 r̄ξ,i r̄γ,i ,i N1 ,R̄x,i : rx,i CzT (r̄ζ,i r̄γ,i ) DT (r̄υ,i r̄ω,i ),i Ñ0 ,s̃0 λs̄ (1 λ)1mI ,ṽ 0 λv̄,κ̃0 (w̃0 )T s̃0 /mI ,where λ [0, 1] is a tuning parameter. When λ 0 the initialpoint becomes the standard starting point, and for λ 1 theinitial point becomes the candidate point.In MPC applications, the optimal control problem is solvedin a receding horizon manner. A good choice of the candidatepoint at time k can therefore be constructed using the solutionfrom the previous time step. As an example consider thesolution of (3) and (4) at time step k 0, for N 3 Tx̂ Tρ Tu Tx̂ Tρ Tu Tx̂ Tρ Tt : u T.011122233The following candidate point is then used at time step k 1 Tx̂ Tρ Tu Tx̂ Tρ Tu Tx̂ Tρ Tt̄ : u T.122233233 r̄γ,N ).Similarly, we left-shift the slack variables, s, and the dualvariables, v and w, to construct s̄, v̄ and w̄.Finally, r̄δ,i : Σ 1δ,i rδ,i . In a similar way we have definedr̄ω,i , r̄η,i , r̄υ,i , r̄ζ,i , r̄ξ,i and r̄γ,i . Solving (10e) for ρi andsubstituting into the remaining equations of (10) yieldsV. P OWER M ANAGEMENT C ASE S TUDYR̄x,N : rx,N CzT (r̄ζ,NB T ṽ0 U0 u0 Ru,0B T ṽi Ui ui Gi x̂i Ru,i ,i Ñ0 x̂1 B u0 Rv,0 x̂i 1 A x̂i B ui Rv,i , ṽi 1 Xi x̂i GTi uiT A ṽi Rx,i ,i Ñ0i Ñ0 ṽN 1 XN x̂N Rx,Nwhere Xi : X̄i Mi Wi 1 MiT and Rx,i : R̄x,i Mi Wi 1 Rw,i . As an example let N 3. In this case, theequations above may be arranged as U0 B T u0Ru,0 ṽ0 Rv,0 B I x̂1 Rx,1 I X1 GT1 AT u1 Ru,1 G1 U1 B T ṽ1 Rv,1 (11) A B I x̂2 Rx,2 I X2 GT2 AT u2 Ru,2 G2 U2 B T A B I ṽ2 Rv,2 x̂3Rx,3 I X3In this section, we compare LPempc against IPMs fromthe following software packages: Gurobi, SeDuMi, MOSEK,LIPSOL, GLPK. These state-of-the-art IPMs are mainly written in low-level language such as FORTRAN and C, and relyon sparse linear algebra that are specifically tailored to thesolution of large-scale sparse linear and conic programs. Wealso include the simplex method provided by CPLEX in ourcomparisons, as well as FORCES that is an IPM based onautomatic code generation [16]. All solvers are called fromMATLAB using MEX interfaces. We have performed oursimulations using an Intel(R) Core(TM) i5-2520M CPU @2.50GHz with 4 GB RAM running a 64-bit Ubuntu 12.04.1LTS operating system.The test system is a system of m generic power generatingunits in the form introduced in [17]. For i 1, 2, . . . , mYi (s) 1(Ui (s) Di (s)) Ei (s).(τi s 1)3(12)Di (s) is the process noise, Ei (s) is the measurement noise,Ui (s) is the power set-point and Yi (s) is the power production.The total productionthe m power generating units is thePm from1sum Z(s) i 1 (τi s 1)3 (Ui (s) Di (s)) . We convert the
SUBMITTED FOR IEEE TRANSACTIONS ON AUTOMATIC CONTROL, AUGUST 20, 2014530Power Plant 1Power Plant 2τipkukuk uk uk903010020000200150-20-402040#iterationsTABLE IC ASE STUDY PARAMETERS2010Plant 1Plant 2z: Totalzmin / mpc (cold started)5001000time [sec]LPempc (warm started)15000.51 0.76 0.883 0.943 0.972 0.986 0.993 0.997 0.998 0.999λ203000030#iterationsproduction [MW]00.51248σ163264128256Fig. 2. Number of iterations as a function of the tuning parameter λ andthe noise parameter σ. Each box-plot has been generated based on an entireclosed-loop simulation. In the top plot we have fixed σ 1, and in the bottomplot λ 0.99.2000210Fig. 1. Closed-loop simulation of a power system controlled by economicMPC. Warm-start (λ 0.99) yields a significant reduction in the number SOLFORCESGLPK010 110110A. SimulationsAn example with two power generating units is considered;a cheap/slow unit, and an expensive/fast unit. This conflictbetween response time and operating costs represents a common situation in the power industry where large thermal powerplants often produce a majority of the electricity, while theuse of units with faster dynamics such as diesel generatorsand gas turbines are limited to critical peak periods. Thecontroller objective is to coordinate the most cost-efficientpower production, respecting capacity constraints and a timevarying electricity demand. It is assumed that full informationabout the initial state is given x0 (0, 0), and that the penaltyof violating the output constraints is qk 104 for all timesteps. Table I lists the system and controller parameters. Weset the prediction horizon to N 80 time steps. It has beenverified by simulation that the controller is stable for this valueof N .Fig. 1 shows a closed-loop simulation with σ 1. Thefigure illustrates the power production of each power generating unit. The cheap unit produces 97% of the energy, whilethe expensive unit is activated only to compensate for thepower imbalances otherwise caused by the slow unit. Fig. 1also shows that warm-start yields a significant reduction in thenumber of iterations. On average the number of iterations isreduced by approximately 37%.horizon length210cpu time [sec]transfer function model into the state space form (1) usinga sampling time of Ts 5 seconds. In the resulting modelstructure, uk Rnu is the nu power set-points, yk Rnyis the measured power production, and zk Rnz is the totalpower production. Note that nu ny m and nz 1. It isassumed that dk N (0, σI) and ek N (0, σI).cpu time [sec]110110010 1101102nr. of units10(a) Increasing N and fixed m 32. (b) Increasing m and fixed N 32.Fig. 3. CPU-time for the different solvers as a function of the horizon Nlength and the number of power generators m.Fig. 2 shows a number of box-plots used to select the warmstart parameter λ. The case λ 0 corresponds to a cold-start.For all values of λ, warm-start reduces the average numberof iterations. We have chosen λ 0.99 for our controller.This value of λ yields an initial point which is both closeto the candidate point and lies well inside the interior of thenon-negative orthant, (w̃, s̃, κ̃, τ̃ ) 0. Fig. 2 shows that forλ 0.99, the number of iterations is reduced even when thevariance of the process and measurement noise is increasedsignificantly.Fig. 3 plots the computation time for solving the LP (2) asa function of the number of power generating units, m, andthe length of the horizon, N . The figure shows that LPmpcis faster than all other solvers with a significant margin. Inaddition to the algorithms included in Fig. 3, the problem (2)was solved using the code generation based IPM CVXGEN[18]. For problems larger than m 4 and N 12, codegeneration in CVXGEN fails due to the problem size. Therefore,we have not included results for CVXGEN in our benchmark.In general, code generation based solvers such as CVXGEN and
SUBMITTED FOR IEEE TRANSACTIONS ON AUTOMATIC CONTROL, AUGUST 20, 20146ACKNOWLEDGEMENTThis work was funded in part by 1) the Danish Ministryof Higher Education and Science in the industrial PhD project”Stochastic MPC with Applications in Smart Energy Systems”(11-117435); 2) the Southern Denmark Growth Forum and theEuropean Regional Development Fund in the project ”Smart& Cool”; and 3) The Danish Council for Strategic Research inthe project ”CITIES - Centre for IT-Intelligent Energy Systemsin Cities” (1305-00027B).1cpu time [sec]10R EFERENCES0100200LPempc (c)400600LPempc (w) Gurobi MOSEK800100012001400time [sec]SeDuMi CPLEX16001800Fig. 4. CPU-time for solving (2) with 15 power generating units and aprediction horizon of 200 time steps. We observe that warm-starting reducesthe average number of iterations by approximately 40%.FORCES are most competitive for small-dimensional problems[16]. Fig. 4 shows CPU-timings for a closed-loop simulationwith 15 power generating units and a prediction horizon ofN 200 time steps. Only the most competitive solvers areincluded. In this simulation LPempc is up to an order ofmagnitude faster than CPLEX, Gurobi, SeDuMi and MOSEK,depending on the problem data. On average, LPempc isapproximately 5 times faster than Gurobi, 6 times faster thanMOSEK, 19 times faster than SeDuMi, and 22 times faster thanCPLEX.VI. C ONCLUSIONSIn this paper, we have developed a computationally efficientIPM for economic MPC of linear systems with a linear objective function and linear constraints. The algorithm combinesthe homogeneous and self-dual model, and a Riccati iterationprocedure specifically tailored to MPC. This is a significantcontribution since existing Riccati iteration procedures forMPC are not directly applicable to the homogeneous and selfdual model that has become widely adopted by state-of-theart IPMs for linear programming. We have also implementedand tested a recent warm-start strategy for homogeneous andself-dual IPMs that has not previously been used in MPCapplications. Our simulations show that this strategy reducesthe average number of iterations by 35-40%, and that aMATLAB and C implementation of the proposed algorithm,LPempc, is significantly faster than several state-of-the-artIPMs, as well as automatic code generation based IPMs forMPC. In a conceptual power management case study, LPempcis up to an order of magnitude faster than CPLEX, Gurobi,SeDuMi and MOSEK. This is important, since the computingtime of solving the optimal control problem is critical inMPC applications. The simulation results also show that thedifference in computing time becomes larger with the problemsize as LPempc scales in a favorable way.[1] S. Mehrotra, “On the Implementation of a Primal-Dual Interior PointMethod,” SIAM Journal on Optimization, vol. 2, no. 4, pp. 575–601,1992.[2] J. Nocedal and S. Wright, Numerical Optimization, ser. Springer Seriesin Operations Research and Financial Engineering. Springer, 2006.[3] A. Skajaa, E. D. Andersen, and Y. Ye, “Warmstarting the homogeneousand self-dual interior point method for linear and conic quadraticproblems,” Mathematical Programming Computation, vol. 5, no. 1, pp.1–25, 2013.[4] C. V. Rao, S. J. Wright, and J. B. Rawlings, “Application of InteriorPoint Methods to Model Predictive Control,” Journal of OptimizationTheory and Applications, vol. 99, no. 3, pp. 723–757, 1998.[5] J. B. Jørgensen, G. Frison, N. F. Gade-Nielsen, and B. Dammann,“Numerical Methods for Solution of the Extended Linear QuadraticControl Problem,” in 4th IFAC Nonlinear Model Predictive ControlConference, 2012, pp. 187–193.[6] J. B. Jørgensen, J. B. Rawlings, and S. B. Jørgensen, “NumericalMethods for Large-Scale Moving Horizon Estimation and Control,” inInternational Symposium on Dynamics and Control Process Systems(DYCOPS), vol. 7, 2004.[7] L. Vandenberghe, S. Boyd, and M. Nouralishahi, “Robust linear programming and optimal control,” in Proceedings of the 15th IFAC WorldCongress, vol. 15, 2002, pp. 271–276.[8] E. D. Andersen, J. Gondzio, C. Meszaros, and X. Xu, “Implementationof Interior Point Methods for Large Scale Linear Programming,” Ecoledes Hautes Etudes Commerciales, Universite de Geneve, Papers 96.3,1996.[9] X. Xu, P.-F. Hung, and Y. Ye, “A simplified homogeneous and selfdual linear programming algorithm and its implementation,” Annals ofOperations Research, vol. 62, no. 1, pp. 151–171, 1996.[10] Y. Ye, M. J. Todd, and S. Mizuno, “An O( nL)-Iteration Homogeneous and Self-Dual Linear Programming Algorithm,” Mathematics ofOperations Research, vol. 19, no. 1, pp. 53–67, 1994.[11] E. D. Andersen, C. Roos, and T. Terlaky, “On implementing a primaldual interior-point method for conic quadratic optimization,” Mathematical Programming, vol. 95, no. 2, pp. 249–277, 2003.[12] V. M. Zavala and L. T. Biegler, “Nonlinear Model Predictive Control,”in Control, Nonlinear Programming Strategies for State Estimationand Model Predictive, ser. Lecture Notes in Control and InformationSciences, L. Magni, D. M. Raimondo, and F. Allgöwer, Eds. SpringerBerlin Heidelberg, 2009, vol. 384, pp. 419–432.[13] M. Diehl, R. Amrit, and J. B. Rawlings, “A Lyapunov Function forEconomic Optimizing Model Predictive Control,” IEEE Transactionson Automatic Control, vol. 56, no. 3, pp. 703–707, 2011.[14] L. Grüne, “Economic receding horizon control without terminal constraints,” Automatica, vol. 49, no. 3, pp. 725–734, 2013.[15] L. E. Sokoler, G. Frison, A. Skajaa, R. Halvgaard, and J. B.Jørgensen, “A Homogeneous and Self-Dual Interior-Point LinearProgramming Algorithm for Economic Model Predictive Control,”Department of Applied Mathematics and Computer Science, TechnicalUniversity of Denmark, Tech. Rep., 2014. [Online]. Available:http://www2.imm.dtu.dk/ jbjo/DTU MPC 2013 25.html[16] A. Domahidi, A. Zgraggen, M. N. Zeilinger, M. Morari, and C. N. Jones,“Efficient Interior Point Methods for Multistage Problems Arising inReceding Horizon Control,” in 51st IEEE Conference on Decision andControl (CDC), 2012, pp. 668–674.[17] K. Edlund, T. Mølbak, and J. D. Bendtsen, “Simple models for modelbased portfolio load balancing controller synthesis,” in 6th IFAC Symposium on Power Plants and Power Systems Control, 2009, pp. 173–178.[18] J. Mattingley and S. Boyd, “CVXGEN: a code generator for embeddedconvex optimization,” Optimization and Engineering, vol. 13, no. 1, pp.1–27, 2012.
the homogeneous and self-dual model can be warm-started efficiently [3].We refer to [8]-[10] for proofs and details. Introduce a new set of optimization variables ( t; v; w; s ), and the additional scalar variables ( ; ). Then the self-dual and homogeneous problem associated with (3) and (4) may be stated as the linear feasibility problem
2 Homogeneous Self-Dual Model Convex optimization has a very strong duality theory that connects the primal problem (1) to its dual, see e.g. [17]. Because of the strong relations between these two problems, modern interior point methods solve simultaneously the two problems making use of information from one to help progress in the other .
data where (1) is denoted dual infeasible if the dual problem max bTy s.t. ATy s c; s 0 (2) corresponding to (1) is infeasible. The vector sis the so-called dual slacks. 2 The homogenous and self dual model However, most methods for solving (1) assume that the problem has an optimal solution. This is in particular true for interior-point methods.
ENUMERATING BASES OF SELF-DUAL MATROIDS Proposition 1.1. Let k be an odd positive integer. If X is an antipodally self-dual cell complex which contains an acyclic, self-dual spanning tree T0, then X gives rise to an involutively self-dual matroid. We use Proposition 1.1, Theorem 1.2 and Lemma 3.2 to obtain the following result. Theorem 1.6.
Avionics: Honeywell 6-Tube EFIS. Enrolled on Honeywell HAPP. Dual Honeywell SPZ-8400 Digital Auto Pilot Dual Honeywell FZ-820 Flight Guidance Computers Dual Honeywell NZ-2010 LR NAV/FMS w/ 6.1 Software Dual Collins VHF-422B COMM (8.33 Spacing) Dual Collins VIR-432 NAV w/ FM Immunity Dual Collins ADF-462 Dual Collins DME-442
1.2.7 Dual monitor set-up screen*9) This page is concerning the dual monitor usage. Use the dual monitor function Check here when you wish to use dual monitor function. OS management dual monitor Select whether it is OS management dual monitor or video card management dual monitor. It is
T. KODA KODAI MATH. J. 10 (1987), 335—342 SELF-DUAL AND ANTI-SELF-DUAL HERMITIAN SURFACES BY TAKASHI KODA 1. Introduction. Let (M, g) be a 4-dimensional oriented Riemannian manifold. The star operator * defined on the space of 2-forms A2M satisfies * *—id. So A2M splits into two eigenspaces as Λ2M Λ2 MQ)Λ2-M, where Λ\M and Λ2-M are the eigenspaces corresponding to eigenvalues 1 and .
self-respect, self-acceptance, self-control, self-doubt, self-deception, self-confidence, self-trust, bargaining with oneself, being one's own worst enemy, and self-denial, for example, are thought to be deeply human possibilities, yet there is no clear agreement about who or what forms the terms between which these relations hold.
Annual Men’s Day Sunday, October 26, 2008 Order of Service Devotion Men’s Day Praise Team Responsive Reading Psalm 1 Congregational Singing This Little Light of Mine Selection Men’s Day Choir Altar Prayer Dance Selection Creations of God Announcements Bro. Robert Shields Pastoral Observations Message Rev. Jeffery A. Lang, Senior Pastor Southside Church & International Ministries Jackson .
