Skip to main content
Log in

A parallel quadratic programming method for dynamic optimization problems

  • Full Length Paper
  • Published:
Mathematical Programming Computation Aims and scope Submit manuscript

Abstract

Quadratic programming problems (QPs) that arise from dynamic optimization problems typically exhibit a very particular structure. We address the ubiquitous case where these QPs are strictly convex and propose a dual Newton strategy that exploits the block-bandedness similarly to an interior-point method. Still, the proposed method features warmstarting capabilities of active-set methods. We give details for an efficient implementation, including tailored numerical linear algebra, step size computation, parallelization, and infeasibility handling. We prove convergence of the algorithm for the considered problem class. A numerical study based on the open-source implementation qpDUNES shows that the algorithm outperforms both well-established general purpose QP solvers as well as state-of-the-art tailored control QP solvers significantly on the considered benchmark problems.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Price excludes VAT (USA)
Tax calculation will be finalised during checkout.

Instant access to the full article PDF.

Fig. 1
Fig. 2
Fig. 3
Fig. 4
Fig. 5
Fig. 6

Similar content being viewed by others

Notes

  1. We make use of the symmetry and positive semidefiniteness of \(\mathcal {M}\, {\tilde{\mathcal {M}}}^{-1} = \mathcal {I}- \delta \, {\tilde{\mathcal {M}}}^{-1}\).

References

  1. qpDUNES Website. http://mathopt.de/qpdunes 2012–2014

  2. Andersson, J.: A general-purpose software framework for dynamic optimization. PhD thesis, Arenberg Doctoral School, KU Leuven, Department of Electrical Engineering (ESAT/SCD) and Optimization in Engineering Center, Kasteelpark Arenberg 10, 3001-Heverlee, Belgium, Oct (2013)

  3. Andersson, J.A.E., Frasch, J.V., Vukov, M., Diehl, M.: A condensing algorithm for nonlinear MPC with a quadratic runtime in horizon length. Automatica (2013, under review)

  4. Bemporad, A., Patrinos, P.: Simple and certifiable quadratic programming algorithms for embedded linear model predictive control. In: Proceedings of the 4th IFAC Nonlinear Model Predictive Control Conference, pp 14–20, (2012)

  5. Berkelaar, A.B., Roos, K., Terlaky, T.: Recent advances in sensitivity analysis and parametric programming. In: The Optimal Set and Optimal Partition Approach to Linear and Quadratic Programming. Kluwer Publishers, Dordrecht (1997)

  6. Bertsekas, D.P., Tsitsiklis, J.N.: Parallel and Distributed Computation: Numerical Methods. Prentice Hall, New Jersey (1989)

    MATH  Google Scholar 

  7. Best M.J.: Applied mathematics and parallel computing. In: An Algorithm for the Solution of the Parametric Quadratic Programming Problem, pp. 57–76. Physica-Verlag, Heidelberg (1996)

  8. Biegler, L.T.: Advances in nonlinear programming concepts for process control. J. Process Control 8(5), 301–311 (1998)

    Article  Google Scholar 

  9. Bock, H.G.: Recent advances in parameter identification techniques for ODE. In: Deuflhard, P., Hairer, E. (eds.) Numerical Treatment of Inverse Problems in Differential and Integral Equations. Birkhäuser, Boston (1983)

    Google Scholar 

  10. Dai, Y.-H., Fletcher, R.: New algorithms for singly linearly constrained quadratic programs subject to low and upper bounds. Math. Program. 106(3), 403–421 (2006)

    Article  MATH  MathSciNet  Google Scholar 

  11. Danskin, J.M.: The Theory of Max-Min and Its Application to Weapons Allocation Problems. Springer-Verlag, Berlin, New York (1967)

    Book  MATH  Google Scholar 

  12. Diehl, M., Bock, H.G., Schlöder, J.P., Findeisen, R., Nagy, Z., Allgöwer, F.: Real-time optimization and nonlinear model predictive control of processes governed by differential-algebraic equations. J. Process Control 12(4), 577–585 (2002)

    Article  Google Scholar 

  13. Domahidi, A., Zgraggen, A., Zeilinger, M.N., Morari, M., Jones, C.N.: Efficient interior point methods for multistage problems arising in receding horizon control. In: IEEE Conference on Decision and Control (CDC), pp. 668–674. Maui (2012)

  14. Ferreau, H.J., Bock, H.G., Diehl, M.: An online active set strategy to overcome the limitations of explicit MPC. Int. J. Robust Nonlinear Control 18(8), 816–830 (2008)

    Article  MATH  MathSciNet  Google Scholar 

  15. Ferreau, H.J., Kirches, C., Potschka, A., Bock, H.G., Diehl, M.: qpOASES: a parametric active-set algorithm for quadratic programming. Math. Program. Comput. (2013)

  16. Ferreau, H.J., Kozma, A., Diehl, M.: A parallel active-set strategy to solve sparse parametric quadratic programs arising in MPC. In: Proceedings of the 4th IFAC Nonlinear Model Predictive Control Conference, pp. 74–79 (2012)

  17. Fiacco, A.V.: Introduction to Sensitivity and Stability Analysis in Nonlinear Programming. Academic Press, New York (1983)

    MATH  Google Scholar 

  18. Fletcher, R.: Practical Methods of Optimization, 2nd edn. Wiley, Chichester (1987)

    MATH  Google Scholar 

  19. Frasch, J.V., Vukov, M., Ferreau, H.J., Diehl, M.: A new quadratic programming strategy for efficient sparsity exploitation in SQP-based nonlinear MPC and MHE. In: Proceedings of the 19th IFAC World Congress (2014)

  20. Frison G., Jorgensen J.: A fast condensing method for solution of linear-quadratic control problems. In: Proceedings of the 52nd IEEE Conference on Decision and Control (2013)

  21. Gerdts, M., Kunkel, M.: A nonsmooth Newton’s method for discretized optimal control problems with state and control constraints. J. Indus. Manage. Optim. 4, 247–270 (2008)

    Article  MATH  MathSciNet  Google Scholar 

  22. Gill, P.E., Murray, W., Wright, M.H.: Practical Optimization. Academic Press, London (1981)

    Google Scholar 

  23. Gill, P.E., Murray, W.: Numerically stable methods for quadratic programming. Math. Program. 14, 349–372 (1978)

    Article  MATH  MathSciNet  Google Scholar 

  24. Houska, B., Ferreau, H.J., Diehl, M.: An auto-generated real-time iteration algorithm for nonlinear MPC in the microsecond range. Automatica 47(10), 2279–2285 (2011)

    Article  MATH  MathSciNet  Google Scholar 

  25. IBM Corp. IBM ILOG CPLEX V12.1, User’s Manual for CPLEX (2009)

  26. Kirches, C., Bock, H.G., Schlöder, J.P., Sager, S.: Block structured quadratic programming for the direct multiple shooting method for optimal control. Optim. Methods Softw. 26, 239–257 (2010)

    Article  Google Scholar 

  27. Kirches, C., Wirsching, L., Bock, H.G., Schlöder, J.P.: Efficient direct multiple shooting for nonlinear model predictive control on long horizons. J. Process Control 22(3), 540–550 (2012)

    Article  MATH  Google Scholar 

  28. Kühl, P., Diehl, M., Kraus, T., Schlöder, J.P., Bock, H.G.: A real-time algorithm for moving horizon state and parameter estimation. Comput. Chem. Eng. 35(1), 71–83 (2011)

    Article  Google Scholar 

  29. Do Duc, L.: Parallelisierung von Innere-Punkte-Verfahren mittels Cyclic Reduction. Bachelor’s thesis, OVGU Magdeburg (2014)

  30. Leineweber, D.B.: Efficient Reduced SQP Methods for the Optimization of Chemical Processes Described by Large Sparse DAE Models. Fortschritt-Berichte VDI Reihe 3, Verfahrenstechnik, vol. 613. VDI Verlag, Düsseldorf (1999)

    Google Scholar 

  31. Li, W., Swetits, J.: A new algorithm for solving strictly convex quadratic programs. SIAM J. Optim. 7(3), 595–619 (1997)

    Article  MATH  MathSciNet  Google Scholar 

  32. Mattingley, J., Boyd, S.: Convex Optimization in Signal Processing and Communications, chapter Automatic Code Generation for Real-Time Convex Optimization. Cambridge University Press, Cambridge (2009)

    Google Scholar 

  33. Mehrotra, S.: On the implementation of a primal-dual interior point method. SIAM J. Optim. 2(4), 575–601 (1992)

    Article  MATH  MathSciNet  Google Scholar 

  34. Nesterov, Y.: Introductory Lectures on Convex Optimization: A Basic Course, vol. 87 of Applied Optimization. Kluwer Academic Publishers, Amsterdam (2004)

    Book  Google Scholar 

  35. Nocedal, J., Wright, S.J.: Numerical Optimization. Springer Series in Operations Research and Financial Engineering, 2nd edn. Springer, New York (2006)

    MATH  Google Scholar 

  36. Qi, L., Sun, J.: A nonsmooth version of Newton’s method. Math. Program. 58, 353–367 (1993)

    Article  MATH  MathSciNet  Google Scholar 

  37. Rao, C.V., Wright, S.J., Rawlings, J.B.: Application of interior-point methods to model predictive control. J. Optim. Theory Appl. 99, 723–757 (1998)

    Article  MATH  MathSciNet  Google Scholar 

  38. Rawlings, J.B., Mayne, D.Q.: Model Predictive Control: Theory and Design. Nob Hill, San Francisco (2009)

    Google Scholar 

  39. Richter, S., Jones, C.N., Morari, M.: Real-time input-constrained MPC using fast gradient methods. In: Proceedings of the IEEE Conference on Decision and Control, Shanghai, China (2009)

  40. Vukov, M., Domahidi, A., Ferreau, H.J., Morari, M., Diehl, M.: Auto-generated algorithms for nonlinear model predicitive control on long and on short horizons. In: Proceedings of the 52nd Conference on Decision and Control (CDC) (2013)

  41. Wang, X., Stoev, J., Pinte, G., Swevers, J.: Energy optimal point-to-point motion using model predictive control. In: Proceedings ASME 2012 5th Annual Dynamic Systems and Control Conference (2012)

  42. Wang, Y., Boyd, S.: Fast model predictive control using online optimization. IEEE Trans. Control Syst. Technol. 18(2), 267–278 (2010)

    Article  Google Scholar 

  43. Wirsching, L., Bock, H.G., Diehl, M.: Fast NMPC of a chain of masses connected by springs. In: Proceedings of the IEEE International Conference on Control Applications, Munich, pp. 591–596 (2006)

  44. Wright, S.: Partitioned dynamic programming for optimal control. SIAM J. Optim. 1(4), 620–642 (1991)

    Article  MATH  MathSciNet  Google Scholar 

  45. Zafiriou, E.: Robust model predictive control of processes with hard constraints. Comput. Chem. Eng. 14(4–5), 359–371 (1990)

    Article  Google Scholar 

Download references

Acknowledgments

This research was supported by Research Council KUL: PFV/10/002 Optimization in Engineering Center OPTEC, GOA/10/09 MaNet and GOA/10/11 Global real-time optimal control of autonomous robots and mechatronic systems. Flemish Government: IOF/KP/SCORES4CHEM, FWO: PhD/postdoc grants and projects: G.0320.08 (convex MPC), G.0377.09 (Mechatronics MPC); IWT: PhD Grants, projects: SBO LeCoPro; Belgian Federal Science Policy Office: IUAP P7 (DYSCO, Dynamical systems, control and optimization, 2012–2017); EU: FP7- EMBOCON (ICT-248940), FP7-SADCO ( MC ITN-264735), ERC ST HIGHWIND (259 166), Eurostars SMART, ACCM.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Janick V. Frasch.

Appendix

Appendix

Proof

(Theorem 6) The proof is done by calculation. We start from the the Cholesky recursion property, apply the assumed relation between Cholesky and Riccati iterates, and transform the expression into the form of the Riccati recursion. We have

$$\begin{aligned} X_{k-1}&= W_{k-1} - U_{k-1} X_k^{-1} U_{k-1}^\top \\ \Leftrightarrow \quad P^{-1} + {\textit{CH}}^{-1} C^\top&= {\textit{CH}}^{-1} C^\top + {\textit{EH}}^{-1} E^\top \\&\quad - {\textit{EH}}^{-1} C^\top \left( P^{-1} + {\textit{CH}}^{-1} C^\top \right) ^{-1} {\textit{CH}}^{-1} E^\top \, , \end{aligned}$$

and therefore

$$\begin{aligned} P^{-1}&= {\textit{EH}}^{-1} E^\top - \textit{EH}^{-1} C^\top \left( P^{-1} + {\textit{CH}}^{-1} C^\top \right) ^{-1} {\textit{CH}}^{-1} E^\top . \end{aligned}$$
(38)

Using the Schur complement \(\bar{Q} = Q - S^\top R^{-1} S\) it is well known from elementary linear algebra that the inverse \(H^{-1}\) can be expressed by

$$\begin{aligned} H^{-1} = \begin{bmatrix} Q&S^\top \\ S&R \end{bmatrix}^{-1} = \begin{bmatrix} \bar{Q}^{-1}&- \bar{Q}^{-1} S^\top R^{-1}\\ - R^{-1} S \bar{Q}^{-1} ~&~ R^{-1} + R^{-1} S \bar{Q}^{-1} S^\top R^{-1} \end{bmatrix}. \end{aligned}$$

Using this, \(C = \begin{bmatrix}\, A&B \,\end{bmatrix}\), and the special structure of \(E = \begin{bmatrix}\, I&0 \,\end{bmatrix}\), we first see that the identities

$$\begin{aligned} {\textit{EH}}^{-1} E^\top&= \bar{Q}^{-1}, \\ {\textit{CH}}^{-1} E^\top&= \left( A - {\textit{BR}}^{-1} S \right) \bar{Q}^{-1} =: \bar{C} \bar{Q}^{-1} \end{aligned}$$

and

$$\begin{aligned} {\textit{CH}}^{-1} C^\top&= \left( A - {\textit{BR}}^{-1} S \right) \bar{Q}^{-1} \left( A^\top - S^\top R^{-1} B^\top \right) + {\textit{BR}}^{-1} B^\top \\&= \bar{C} \bar{Q}^{-1} \bar{C}^\top + {\textit{BR}}^{-1} B^\top \end{aligned}$$

hold. Thus, (38) can be written as

$$\begin{aligned} P^{-1}&= \bar{Q}^{-1} - \bar{Q}^{-1} \bar{C}^\top \left( P^{-1} + B R^{-1} B^\top + \bar{C} \bar{Q}^{-1} \bar{C}^\top \right) ^{-1} \bar{C} \bar{Q}^{-1} \, . \end{aligned}$$

Applying the Woodbury matrix identity with \(Y := P^{-1} + {\textit{BR}}^{-1} B^\top \) we can express the right-hand-side term by

$$\begin{aligned} P^{-1} = \left( \bar{Q} + \bar{C}^\top Y^{-1} \bar{C} \right) ^{-1} \,. \end{aligned}$$

Thus

$$\begin{aligned} P&= \bar{Q} + \bar{C}^\top Y^{-1} \bar{C} \\&= Q - S^\top R^{-1} S + \left( A^\top - S^\top R^{-1} B^\top \right) \left( P^{-1} + {\textit{BR}}^{-1} B^\top \right) ^{-1} \left( A - {\textit{BR}}^{-1} S \right) . \end{aligned}$$

Note that the inverse matrices of \(Q,R,P\) and \(Y\) exist (and are real-valued), since we assume that QP (1) is convex. Applying the Woodbury identity once again (however, in opposite direction) on \(\left( P^{-1} + {\textit{BR}}^{-1} B^\top \right) ^{-1}\), and introducing \(\bar{R} := \left( R + B^\top {\textit{PB}} \right) \), we get

$$\begin{aligned} P&= Q - S^\top R^{-1} S \nonumber \\&\quad + \left( A^\top - S^\top R^{-1} B^\top \right) \left( P - {\textit{PB}} \left( R + B^\top {\textit{PB}} \right) ^{-1} B^\top P \right) \left( A - {\textit{BR}}^{-1} S \right) \nonumber \\&= Q - S^\top R^{-1} S \nonumber \\&\quad + A^\top \left( P - {\textit{PB}} \bar{R}^{-1} B^\top P \right) A \nonumber \\&\quad - S^\top R^{-1} B^\top \left( P - {\textit{PB}} \bar{R}^{-1} B^\top P \right) A \nonumber \\&\quad - A^\top \left( P - {\textit{PB}} \bar{R}^{-1} B^\top P \right) {\textit{BR}}^{-1} S \nonumber \\&\quad + S^\top R^{-1} B^\top \left( P - {\textit{PB}} \bar{R}^{-1} B^\top P \right) {\textit{BR}}^{-1} S. \end{aligned}$$
(39)

Using the identity \(I = \bar{R} \bar{R}^{-1} = \bar{R}^{-1} \bar{R}\), we further have

$$\begin{aligned}&S^\top R^{-1} B^\top \left( P - {\textit{PB}} \bar{R}^{-1} B^\top P \right) {\textit{BR}}^{-1} S - S^\top R^{-1} S \\&\quad = S^\top R^{-1} B^\top P {\textit{BR}}^{-1} S - S^\top R^{-1} B^\top {\textit{PB}} \bar{R}^{-1} B^\top P {\textit{BR}}^{-1} S - S^\top R^{-1} S \\&\quad = S^\top R^{-1} \bar{R} \bar{R}^{-1} B^\top P {\textit{BR}}^{-1} S \\&\qquad - S^\top R^{-1} B^\top {\textit{PB}} \bar{R}^{-1} B^\top P {\textit{BR}}^{-1} S - S^\top \bar{R}^{-1} \bar{R} R^{-1} S \\&\quad = S^\top R^{-1} \left( R + B^\top {\textit{PB}} \right) \bar{R}^{-1} B^\top P {\textit{BR}}^{-1} S \\&\qquad - S^\top R^{-1} B^\top {\textit{PB}} \bar{R}^{-1} B^\top P {\textit{BR}}^{-1} S - S^\top \bar{R}^{-1} \left( R + B^\top {\textit{PB}} \right) R^{-1} S \\&\quad = S^\top R^{-1} R \bar{R}^{-1} B^\top P {\textit{BR}}^{-1} S + S^\top R^{-1} B^\top {\textit{PB}} \bar{R}^{-1} B^\top P {\textit{BR}}^{-1} S \\&\qquad - S^\top R^{-1} B^\top {\textit{PB}} \bar{R}^{-1} B^\top P {\textit{BR}}^{-1} S \\&\qquad - S^\top \bar{R}^{-1} R R^{-1} S - S^\top \bar{R}^{-1} B^\top {\textit{PB}} R^{-1} S \\&\quad = - S^\top \bar{R}^{-1} S \end{aligned}$$

and

$$\begin{aligned}&- A^\top \left( P - {\textit{PB}} \bar{R}^{-1} B^\top P \right) {\textit{BR}}^{-1} S \\&\quad = - A^\top P {\textit{BR}}^{-1} S + A^\top {\textit{PB}} \bar{R}^{-1} B^\top P {\textit{BR}}^{-1} S \\&\quad = - A^\top {\textit{PB}} \bar{R}^{-1} \bar{R} R^{-1} S + A^\top {\textit{PB}} \bar{R}^{-1} B^\top P {\textit{BR}}^{-1} S \\&\quad = - A^\top {\textit{PB}} \bar{R}^{-1} \left( R + B^\top {\textit{PB}} \right) R^{-1} S + A^\top {\textit{PB}} \bar{R}^{-1} B^\top P {\textit{BR}}^{-1} S \\&\quad = - A^\top {\textit{PB}} \bar{R}^{-1} {\textit{RR}}^{-1} S - A^\top {\textit{PB}} \bar{R}^{-1} B^\top P {\textit{BR}}^{-1} S + A^\top {\textit{PB}} \bar{R}^{-1} B^\top P {\textit{BR}}^{-1} S \\&\quad = - A^\top {\textit{PB}} \bar{R}^{-1} S . \end{aligned}$$

Analogously,

$$\begin{aligned} - S^\top R^{-1} B^\top \left( P - {\textit{PB}} \bar{R}^{-1} B^\top P \right) A = - S^\top \bar{R}^{-1} B^\top P A . \end{aligned}$$

Therefore (39) is equivalent to

$$\begin{aligned} P&= Q - S^\top \bar{R}^{-1} S \\&\quad - S^\top \bar{R}^{-1} B^\top P A - A^\top {\textit{PB}} \bar{R}^{-1} S + A^\top \left( P - {\textit{PB}} \bar{R}^{-1} B^\top P \right) A \\&= Q+ A^\top P A - \left( S^\top + A^\top {\textit{PB}} \right) \bar{R}^{-1} \left( S + B^\top P A \right) , \end{aligned}$$

which concludes the proof. \(\square \)

Proof

(Corollary 7 ) The proof follows exactly the lines of the proof to Theorem 6, but keeps the matrix block indices. In particular, one has

$$\begin{aligned} P_{k-1}^{-1}&= E_{k-1} H_{k-1}^{-1} E_{k-1}^\top - E_{k-1} H_{k-1}^{-1} C_{k-1}^\top \left( P_{k}^{-1} + C_{k-1} H_{k-1}^{-1} C_{k-1}^\top \right) ^{-1}\\&\quad \times C_{k-1} H_{k-1}^{-1} E_{k-1}^\top \end{aligned}$$

in place of (38) and transforms it into (19b) using the same matrix identities. \(\square \)

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Frasch, J.V., Sager, S. & Diehl, M. A parallel quadratic programming method for dynamic optimization problems. Math. Prog. Comp. 7, 289–329 (2015). https://doi.org/10.1007/s12532-015-0081-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1007/s12532-015-0081-7

Keywords

Mathematics Subject Classification

Navigation