0018-9286 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAC.2015.2495558, IEEE

Transactions on Automatic Control

SUBMITTED FOR IEEE TRANSACTIONS ON AUTOMATIC CONTROL, SEPTEMBER 30, 2015 1

A Homogeneous and Self-Dual Interior-Point

Linear Programming Algorithm for Economic

Model Predictive Control

Leo Emil Sokoler, Gianluca Frison, Anders Skajaa, Rasmus Halvgaard, John Bagterp Jørgensen

Abstract—We develop an efficient homogeneous and self-dual interior-point method (IPM) for the linear programs arising in economic model predictive control of constrained linear systems with linear objective functions. The algorithm is based on a

Riccati iteration procedure, which is adapted to the linear system of equations solved in homogeneous and self-dual IPMs. Fast convergence is further achieved using a warm-start strategy. We implement the algorithm in MATLAB and C. Its performance is tested using a conceptual power management case study.

Closed loop simulations show that 1) the proposed algorithm is significantly faster than several state-of-the-art IPMs based on sparse linear algebra, and 2) warm-start reduces the average number of iterations by 35-40%.

Index Terms—Optimization algorithms, Linear programming algorithms, Predictive control for linear systems, Riccati iterations, Energy systems

I. INTRODUCTION

In economic model predictive control (EMPC) of linear systems with a linear objective function and linear constraints, the optimal control problem (OCP) can be posed as a linear program (LP). As the optimization problem is solved online, the performance and reliability of the optimization algorithm solving the LP is important. The homogeneous and self-dual model has become widely adopted by state-of-the-art IPMs to solve LPs. This paper presents a homogeneous and self-dual variant of Mehrotra’s predictor-corrector method [1], [2] for

EMPC. The algorithm combines the following methods for computational efficiency: • Riccati Iteration Procedure: The most time consuming numerical operations in the IPM are handled by a Riccati iteration procedure tailored to EMPC. Riccati iteration procedures have not previously been combined with the homogeneous and self-dual model. • Warm-Start: The warm-start strategy of [3] is applied to reduce the number of IPM iterations. This strategy, which is designed for homogeneous and self-dual IPMs, has not been applied to MPC problems in the existing literature.

In [3] warm-start reduces the number of iterations by 3075% based on the NETLIB collection of test problems.

Riccati-based IPMs have been developed for set-point based

MPC with an `2-penalty [4], [5] and with an `1-penalty [6].

The main contributions of this paper are 1) to combine the

The authors are with the Department of Applied Mathematics and Computer

Science, DTU, DK-2800 Kgs. Lyngby, Denmark. L. E. Sokoler and A. Skajaa are also affiliated with DONG Energy, DK-2830 Virum, Denmark (e-mail: {leso, giaf, andsk, rhal, jbjo}@dtu.dk) homogeneous and self dual model with a Riccati iteration procedure, 2) to test the warm-start procedure of [3] in an

MPC framework, and 3) to report performance results for an efficient implementation of the proposed algorithm.

This paper is organized as follows. Section II formulates the

OCP solved in EMPC as a highly structured LP. Section III develops a special purpose homogeneous and self-dual IPM for EMPC. Warm-start is introduced in Section IV. Section V compares a MATLAB and C implementation of the proposed algorithm denoted LPempc to several state-of-the-art IPMs using a simple power management case study1. Concluding remarks are given in Section VI.

II. OPTIMAL CONTROL PROBLEM

We consider linear state space systems in the form xk+1 = Axk +Buk + Edk, dk ∼ N(0, Rd), (1a) yk = Cyxk + ek, ek ∼ N(0, Re), (1b) zk = Czxk, (1c) where x0 ∼ N(xˆ0, P0). (A,B,Cy, Cz, E) are the state space matrices, xk ∈ Rnx is the state vector, uk ∈ Rnu is the input vector, yk ∈ Rny is the measurement vector, zk ∈ Rnz is the output vector, dk is the process noise vector. and ek is the measurement noise vector. We use bold letters to denote stochastic variables. In this paper, the OCP solved at every sampling instant is defined as min. u,xˆ,zˆ,ρ ∑ j∈N0 pTk+juk+j + q

T k+j+1ρk+j+1, (2a) s.t. xˆk+j+1|k = Axˆk+j|k +Buk+j , j ∈ N0, (2b) zˆk+j|k = Czxˆ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) zk+j − ρk+j ≤ zˆk+j|k ≤ zk+j + ρk+j , j ∈ N1, (2f) ρk+j ≥ 0, j ∈ N1, (2g) where Ni := {0 + i, 1 + i, . . . , N − 1 + i}, with N being the length of the prediction horizon. The problem data are the state-space matrices, (A,B,Cz), the filtered estimate, xˆk|k, the input limits, (uk+j , uk+j), the input-rate limits, (∆uk+j ,∆uk+j), the output limits, (zk+j , zk+j), the input 1Software is available via http://www2.imm.dtu.dk/∼jbjo/publications 0018-9286 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAC.2015.2495558, IEEE

Transactions on Automatic Control

SUBMITTED FOR IEEE TRANSACTIONS ON AUTOMATIC CONTROL, SEPTEMBER 30, 2015 2 prices, pk+j , and the price for violating the output constraints, qk+j . We use soft output constraints, defined by (2f) and (2g), to ensure feasibility of the OCP. For compact notation, the optimization variables in (2) are written as u, xˆ, zˆ and ρ, where u = [ uTk u