Publications
Journal Articles
2018

The Voice of OptimizationArXiv eprints, Dec. 2018
We introduce the idea that using optimal classification trees (OCTs) and optimal classification trees withhyperplanes (OCTHs), interpretable machine learning algorithms developed in [BD17, BD18], we are able to obtain insight on the strategy behind the optimal solution in any continuous and mixedinteger convex optimization problem as a function of key parameters that affect the problem. In this way, optimization is not a black box anymore. Instead, we redefine optimization as a multiclass classification problem where the predictor gives insights on the logic behind the optimal solution. In other words, OCTs and OCTHs give optimization a voice. We show on several realistic examples that the accuracy behind our method is in the 90%100% range, while even when the predictions are not correct, the degree of suboptimality or infeasibility is very low. We compare optimal strategy predictions of OCTs and OCTHs and feedforward neural networks (NNs) and conclude that the performance of OCTHs and NNs is comparable. OCTs are somewhat weaker but often competitive. Therefore, our approach provides a novel, reliable and insightful understanding of optimal strategies to solve a broad class of continuous and mixedinteger optimization problems.
@article{stellato2018b, author = {{Bertsimas}, D. and {Stellato}, B.}, title = {The Voice of Optimization}, journal = {ArXiv eprints}, year = {2018}, month = dec, adsnote = {Provided by the SAO/NASA Astrophysics Data System}, adsurl = {https://ui.adsabs.harvard.edu/\#abs/2018arXiv181209991B}, archiveprefix = {arXiv}, eprint = {1812.09991}, keywords = {Mathematics  Optimization and Control}, pdf = {https://arxiv.org/pdf/1812.09991.pdf}, primaryclass = {math.OC}, url = {https://arxiv.org/pdf/1812.09991} }

OSQP: An Operator Splitting Solver for Quadratic ProgramsArXiv eprints, Jan. 2018
We present a general purpose solver for quadratic programs based on the alternating direction method of multipliers, employing a novel operator splitting technique that requires the solution of a quasidefinite linear system with the same coefficient matrix in each iteration. Our algorithm is very robust, placing no requirements on the problem data such as positive definiteness of the objective function or linear independence of the constraint functions. It is divisionfree once an initial matrix factorization is carried out, making it suitable for realtime applications in embedded systems. In addition, our technique is the first operator splitting method for quadratic programs able to reliably detect primal and dual infeasible problems from the algorithm iterates. The method also supports factorization caching and warm starting, making it particularly efficient when solving parametrized problems arising in finance, control, and machine learning. Our opensource C implementation OSQP has a small footprint, is libraryfree, and has been extensively tested on many problem instances from a wide variety of application areas. It is typically ten times faster than competing interior point methods, and sometimes much more when factorization caching or warm start is used.
@article{stellato2018a, author = {{Stellato}, B. and {Banjac}, G. and {Goulart}, P. and {Bemporad}, A. and {Boyd}, S.}, title = {{OSQP: An Operator Splitting Solver for Quadratic Programs}}, journal = {ArXiv eprints}, year = {2018}, month = jan, adsnote = {Provided by the SAO/NASA Astrophysics Data System}, adsurl = {http://adsabs.harvard.edu/abs/2017arXiv171108013S}, archiveprefix = {arXiv}, code = {https://osqp.org}, eprint = {1711.08013}, keywords = {Mathematics  Optimization and Control}, pdf = {https://arxiv.org/pdf/1711.08013.pdf}, primaryclass = {math.OC}, url = {https://arxiv.org/abs/1711.08013} }
2017

SecondOrder Switching Time Optimization for Switched Dynamical SystemsIEEE Transactions on Automatic Control, Vol. 62, no. 10, pp. 5407–5414, Oct. 2017
Switching time optimization arises in finitehorizon optimal control for switched systems where, given a sequence of continuous dynamics, one minimizes a cost function with respect to the switching times. We propose an efficient method for computing the optimal switching times for switched linear and nonlinear systems. A novel secondorder optimization algorithm is introduced where, at each iteration, the dynamics are linearized over an underlying time grid to compute the cost function, the gradient and the Hessian efficiently. With the proposed method, the most expensive operations at each iteration are shared between the cost function and its derivatives, thereby greatly reducing the computational burden. We implemented the algorithm in the Julia package SwitchTimeOpt allowing the user to easily solve switching time optimization problems. In the case of linear dynamics, many operations can be further simplified and benchmarks show that our approach is able to provide optimal solutions in just a few ms. In the case of nonlinear dynamics, two examples show that our method provides optimal solutions with up to two orders of magnitude time reductions over stateoftheart approaches.
@article{stellato2017a, author = {Stellato, B. and OberBl{\"o}baum, S. and Goulart, P.}, title = {SecondOrder Switching Time Optimization for Switched Dynamical Systems}, journal = {{IEEE} Transactions on Automatic Control}, year = {2017}, volume = {62}, number = {10}, pages = {54075414}, month = oct, issn = {00189286}, code = {https://github.com/oxfordcontrol/SwitchTimeOpt.jl}, doi = {10.1109/TAC.2017.2697681}, groups = {switching time optimization}, pdf = {http://arxiv.org/pdf/1608.08597.pdf}, url = {https://doi.org/10.1109/TAC.2017.2697681} }

Infeasibility detection in the alternating direction method of multipliers for convex optimizationoptimizationonline.org, Jun. 2017
The alternating direction method of multipliers (ADMM) is a powerful operator splitting technique for solving structured optimization problems. For convex optimization problems, it is well known that the iterates generated by ADMM converge to a solution provided that it exists. If a solution does not exist, then the ADMM iterates do not converge. Nevertheless, we show that the ADMM iterates yield conclusive information regarding problem infeasibility for a wide class of convex optimization problems including both quadratic and conic programs. In particular, we show that in the limit the ADMM iterates either satisfy a set of firstorder optimality conditions or produce a certificate of either primal or dual infeasibility. Based on these results, we propose termination criteria for detecting primal and dual infeasibility in ADMM.
@article{banjac2017, author = {Banjac, G. and Goulart, P. and Stellato, B. and Boyd, S.}, title = {Infeasibility detection in the alternating direction method of multipliers for convex optimization}, journal = {optimizationonline.org}, year = {2017}, month = jun, groups = {quadratic programs}, pdf = {http://www.optimizationonline.org/DB_FILE/2017/06/6058.pdf}, url = {http://www.optimizationonline.org/DB_HTML/2017/06/6058.html} }

HighSpeed Finite Control Set Model Predictive Control for Power ElectronicsIEEE Transactions on Power Electronics, Vol. 32, no. 5, pp. 4007–4020, May. 2017First Prize Paper Award IEEE Transactions on Power Electronics
Common approaches for direct model predictive control (MPC) for current reference tracking in power electronics suffer from the high computational complexity encountered when solving integer optimal control problems over long prediction horizons. We propose an efficient alternative method based on approximate dynamic programming, greatly reducing the computational burden and enabling sampling times under \(25\: μs\). Our approach is based on the offline minimization of an infinite horizon value function estimate which is then applied to the tail cost of the MPC problem. This allows us to reduce the controller horizon to a very small number of stages while improving overall controller performance. Our proposed algorithm is implemented on a small size FPGA and validated on a variable speed drive system with a threelevel voltage source converter. Time measurements showed that only \(5.76\: μs\)are required to run our algorithm for horizon N=1 and \(17.27\: μs\)for N=2 while outperforming state of the art approaches with much longer horizons in terms of currents distortion and switching frequency. To the authors’ knowledge, this is the first time direct MPC for current control has been implemented on an FPGA solving the integer optimization problem in realtime and achieving comparable performance to formulations with long prediction horizons.
@article{stellato2017, author = {Stellato, B. and Geyer, T. and Goulart, P.}, title = {HighSpeed Finite Control Set Model Predictive Control for Power Electronics}, journal = {{IEEE} Transactions on Power Electronics}, year = {2017}, volume = {32}, number = {5}, pages = {40074020}, month = may, issn = {08858993}, prize = {First Prize Paper Award IEEE Transactions on Power Electronics}, prize_link = {https://www.ieeepels.org/awards/ieeetransactionsonpowerelectronicsprizepaperaward}, groups = {power electronics}, pdf = {http://arxiv.org/pdf/1510.05578.pdf}, url = {http://dx.doi.org/10.1109/TPEL.2016.2584678} }

Multivariate Chebyshev Inequality with Estimated Mean and VarianceThe American Statistician, Vol. 71, no. 2, pp. 123–127, 2017
A variant of the wellknown Chebyshev inequality for scalar random variables can be formulated in the case where the mean and variance are estimated from samples. In this paper we present a generalization of this result to multiple dimensions where the only requirement is that the samples are independent and identically distributed. Furthermore, we show that as the number of samples tends to infinity our inequality converges to the theoretical multidimensional Chebyshev bound.
@article{stellato2017b, author = {Stellato, B. and Van Parys, B. P.G. and Goulart, P.}, title = {Multivariate Chebyshev Inequality with Estimated Mean and Variance}, journal = {The American Statistician}, year = {2017}, volume = {71}, number = {2}, pages = {123127}, groups = {probability inequalities}, pdf = {http://arxiv.org/pdf/1509.08398.pdf}, url = {http://dx.doi.org/10.1080/00031305.2016.1186559} }
Conference Proceedings
2018

Embedded MixedInteger Quadratic Optimization Using the OSQP SolverEuropean Control Conference (ECC), Jul. 2018
We present a novel branchandbound solver for mixedinteger quadratic programs (MIQPs) that efficiently exploits the firstorder OSQP solver for the quadratic program (QP) subproblems. Our algorithm is very robust, requires no dynamic memory allocation and is divisionfree once an initial factorization is computed. Thus, it suitable for embedded applications with low computing power. Moreover, it does not require any assumption on the problem data such as strict convexity of the objective function. We exploit factorization caching and warmstarting to reduce the computational cost of QP relaxations during branchandbound and over repeated solutions of parametric MIQPs such as those arising in embedded control, portfolio optimization, and machine learning. Numerical examples show that our method, using a simple highlevel Python implementation interfaced with the OSQP solver, is competitive with established commercial solvers.
@inproceedings{stellato2018, author = {Stellato, B. and Naik, V. V. and Bemporad, A. and Goulart, P. and Boyd, S.}, title = {Embedded MixedInteger Quadratic Optimization Using the {OSQP} Solver}, booktitle = {European Control Conference ({ECC})}, year = {2018}, code = {https://github.com/oxfordcontrol/miosqp}, month = jul, pdf = {/assets/downloads/publications/2018/miosqp_ecc.pdf}, groups = {power electronics, integer programs} }
2017

Embedded Code Generation Using the OSQP SolverIEEE Conference on Decision and Control (CDC), Dec. 2017
We introduce a code generation software package that accepts a parametric description of a quadratic program (QP) as input and generates tailored C code that compiles into a fast and reliable optimization solver for the QP that can run on embedded platforms. The generated code is based on OSQP, a novel opensource operator splitting solver for quadratic programming. Our software supports matrix factorization caching and warm starting, and allows updates of the problem parameters during runtime. The generated C code is libraryfree and has a very small compiled footprint. Examples arising in realworld applications show that the generated code outperforms stateoftheart embedded and desktop QP solvers.
@inproceedings{banjac2017a, author = {Banjac, G. and Stellato, B. and Moehle, N. and Goulart, P. and Bemporad, A. and Boyd, S.}, title = {Embedded Code Generation Using the {OSQP} Solver}, booktitle = {{IEEE} Conference on Decision and Control ({CDC})}, year = {2017}, month = dec, pdf = {/assets/downloads/publications/2017/osqp_embedded.pdf}, url = {https://ieeexplore.ieee.org/document/8263928}, groups = {quadratic programs} }
2016

Realtime FPGA Implementation of Direct MPC for Power ElectronicsIEEE Conference on Decision and Control (CDC), pp. 1471–1476, Dec. 2016
Common approaches for direct model predictive control (MPC) for current reference tracking in power electronics suffer from the high computational complexity encountered when solving integer optimal control problems over long prediction horizons. Recently, an alternative method based on approximate dynamic programming showed that it is possible to reduce the computational burden enabling sampling times under \(25\:μs\)by shortening the MPC horizon to a very small number of stages while improving the overall controller performance. In this paper we implemented this new approach on a small size FPGA and validated it on a variable speed drive system with a threelevel voltage source converter. Time measurements showed that only \(5.76\:μs\)are required to run our algorithm for horizon \(N=1\)and \(17.27\:μs\)for \(N=2\)while outperforming state of the art approaches with much longer horizons in terms of currents distortion and switching frequency. To the authors knowledge, this is the first time direct MPC for current control has been implemented on an embedded platform achieving comparable performance to formulations with long prediction horizons.
@inproceedings{stellato2016a, author = {Stellato, B. and Goulart, P.}, title = {Realtime {FPGA} Implementation of Direct {MPC} for Power Electronics}, booktitle = {{IEEE} Conference on Decision and Control ({CDC})}, year = {2016}, pages = {14711476}, month = dec, groups = {power electronics}, pdf = {/assets/downloads/publications/2016/pels_cdc2016.pdf}, url = {https://doi.org/10.1109/CDC.2016.7798474} }

Optimal Control of Switching Times in Switched Linear SystemsIEEE Conference on Decision and Control (CDC), pp. 7228–7233, Dec. 2016
Switching time optimization arises in finitehorizon optimal control for switched systems where, given a sequence of continuous dynamics, we minimize a cost function with respect to the switching times. In this paper we propose an efficient method for computing optimal switching times in switched linear systems. We derive simple expressions for the cost function, the gradient and the Hessian which can be computed efficiently online without performing any integration. With the proposed method, the most expensive computations are decomposed into independent scalar exponentials which can be efficiently computed and parallelized. Simulation results show that our method is able to provide fast convergence and handle efficiently a high number switching times.
@inproceedings{stellato2016b, author = {Stellato, B. and OberBl\"{o}baum, S. and Goulart, P.}, title = {Optimal Control of Switching Times in Switched Linear Systems}, booktitle = {{IEEE} Conference on Decision and Control ({CDC})}, year = {2016}, pages = {72287233}, month = dec, groups = {switching time optimization}, pdf = {/assets/downloads/publications/2016/switchtimeopt2016.pdf}, url = {https://doi.org/10.1109/CDC.2016.7799384} }

HighSpeed Direct Model Predictive Control for Power ElectronicsEuropean Control Conference (ECC), pp. 129–134, Jul. 2016
Common approaches for direct model predictive control (MPC) for current reference tracking in power electronics suffer from the high computational complexity encountered when solving integer optimal control problems over long prediction horizons. We propose an efficient alternative method based on approximate dynamic programming, greatly reducing the computational burden and enabling sampling times under \(25\:μs\). Our approach is based on the offline minimization of an infinite horizon cost function estimate which is then applied to the tail cost of the MPC problem. This allows us to reduce the controller horizon to a very small number of stages improving overall controller performance. Our proposed algorithm is validated on a variable speed drive system with a threelevel voltage source converter.
@inproceedings{stellato2016, author = {Stellato, B. and Goulart, P.}, title = {HighSpeed Direct Model Predictive Control for Power Electronics}, booktitle = {European Control Conference ({ECC})}, year = {2016}, pages = {129134}, month = jul, groups = {power electronics}, pdf = {/assets/downloads/publications/2016/pel_ecc2016.pdf}, url = {http://ieeexplore.ieee.org/document/7810275/} }
Theses
2017

MixedInteger Optimal Control of Fast Dynamical SystemsUniversity of Oxford, 2017
Many applications in engineering, computer science and economics involve mixedinteger optimal control problems. Solving these problems in realtime is a challenging task because of the explosion of integer combinations to evaluate. This thesis focuses on the development of new algorithms for mixedinteger programming with an emphasis on optimal control problems of fast dynamical systems with discrete controls. The first part proposes two reformulations to reduce the computational complexity. The first reformulation avoids integer variables altogether. By considering a sequence of switched dynamics, we analyze the switching time optimization problem. Even though it is a continuous smooth problem, it is nonconvex and the cost function and derivatives are hard to compute. We develop a new efficient method to compute the cost function and its derivatives. Our technique brings up to two orders of magnitude speedups with respect to stateoftheart tools. The second approach reduces the number of integer decisions. In hybrid model predictive control (MPC) the computational complexity grows exponentially with the horizon length. Using approximate dynamic programming (ADP) we reduce the horizon length while maintaining good control performance by approximating the tail cost offline. This approach allows, for the first time, the application of such control techniques to fast dynamical systems with sampling times of only a few microseconds. The second part investigates embedded branchandbound algorithms for mixedinteger quadratic programs (MIQPs). A core component of these methods is the solution of continuous quadratic programs (QPs). We develop OSQP, a new robust and efficient generalpurpose QP solver based on the alternating direction method of multipliers (ADMM) and able, for the first time, to detect infeasible problems. We include OSQP into a custom branchandbound algorithm suitable for embedded systems. Our extension requires only a single matrix factorization and exploits warmstarting, thereby greatly reducing the number of ADMM iterations required. Numerical examples show that our algorithm solves small to medium scale MIQPs more quickly than commercial solvers.
@phdthesis{stellato2017d, author = {Stellato, B.}, title = {MixedInteger Optimal Control of Fast Dynamical Systems}, school = {University of Oxford}, year = {2017}, groups = {theses}, owner = {sidereus}, pdf = {/assets/downloads/publications/stellato_thesis.pdf}, url = {https://ora.ox.ac.uk/objects/uuid:b8a7323ce36e45ecae8d6c9eb4350629}, timestamp = {2017.12.20} }
2014

Datadriven chance constrained optimizationETH Zürich, 2014
Traditional optimization methods for decisonmaking under uncertainty assume perfect model information. In practice, however, such precise knowledge is rarely available. Thus, optimal decisions coming from these approaches can be very sensitive to perturbations and unreliable. Stochastic optimization programs take into account uncertainty but are intractable in general and need to be approximated. Of late, distributionally robust optimization methods have shown to be powerful tools to reformulate stochastic programs in a tractable way. Moreover, the recent advent of cheap sensing devices has caused the explosion of available historical data, usually referred to as “Big Data”. Thus, modern optimization techniques are shifting from traditional methods to datadriven approaches. In this thesis, we derive datadriven tractable reformulations for stochastic optimization programs based on distributionally robust optimization. In the first part of this work we provide our theoretical contributions. New distributionally robust probability bounds are derived and used to reformulate uncertain optimization programs assuming limited information about the uncertainty. Then, we show how this information can be derived from historical data. In the second part of this work, we compare the developed methods to support vector machines in a machine learning setting and to randomized optimization and in a control context.
@mastersthesis{stellato2014, author = {Stellato, B.}, title = {Datadriven chance constrained optimization}, school = {ETH Z\"{u}rich}, year = {2014}, groups = {theses}, pdf = {/assets/downloads/publications/stellato_master_thesis.pdf}, url = {http://dx.doi.org/10.3929/ethza010266857} }