Variable Horizon Model Predictive
Control: Robustness & Optimality
16th April 2012
(Revised 8th June 2012)
Rohan C. Shekhar
Pembroke College
Control Group
Department of Engineering
University of Cambridge
A thesis submitted for the degree of
Doctor of Philosophy

Variable Horizon Model Predictive Control (VH-MPC) is a form of predictive control that
includes the horizon length as a decision variable in the constrained optimisation problem
solved at each iteration. It has been recently applied to completion problems, where the system
state is to be steered to a closed set in finite time. The behaviour of the system once completion
has occurred is not considered part of the control problem.
This thesis is concerned with three aspects of robustness and optimality in VH-MPC completion
problems. In particular, the thesis investigates robustness to well defined but unpredictable
changes in system and controller parameters, robustness to bounded disturbances in the
presence of certain input parameterisations to reduce computational complexity, and optimal
robustness to bounded disturbances using tightened constraints. In the context of linear time
invariant systems, new theoretical contributions and algorithms are developed.
Firstly, changing dynamics, constraints and control objectives are addressed by introducing the
notion of feasible contingencies. A novel algorithm is proposed that introduces extra prediction
variables to ensure that anticipated new control objectives are always feasible, under changed
system parameters. In addition, a modified constraint tightening formulation is introduced
to provide robust completion in the presence of bounded disturbances. Different contingency
scenarios are presented and numerical simulations demonstrate the formulation’s efficacy.
Next, complexity reduction is considered, using a form of input parameterisation known as
move blocking. After introducing a new notation for move blocking, algorithms are presented
for designing a move-blocked VH-MPC controller. Constraints are tightened in a novel way for
robustness, whilst ensuring that guarantees of recursive feasibility and finite-time completion
are preserved. Simulations are used to illustrate the effect of an example blocking scheme on
computation time, closed-loop cost, control inputs and state trajectories.
Attention is now turned towards mitigating the effect of constraint tightening policies on a
VH-MPC controller’s region of attraction. An optimisation problem is formulated to maximise
the volume of an inner approximation to the region of attraction, parameterised in terms of
the tightening policy. Alternative heuristic approaches are also proposed to deal with high
state dimensions. Numerical examples show that the new technique produces substantially
improved regions of attraction in comparison to other proposed approaches, and greatly
reduces the maximum required prediction horizon length for a given application.
Finally, a case study is presented to illustrate the application of the new theory developed in
this thesis to a non-trivial example system. A simplified nonlinear surface excavation machine
and material model is developed for this purpose. The model is stabilised with an inner-loop
controller, following which a VH-MPC controller for autonomous trajectory generation is
designed using a discretised, linearised model of the stabilised system. Realistic simulated
trajectories are obtained from applying the controller to the stabilised system and incorporating
the ideas developed in this thesis.
These ideas improve the applicability and computational tractability of VH-MPC, for both
traditional applications as well as those that go beyond the realm of vehicle manœuvring.
Abstract
Variable Horizon Model Predictive Control: Robustness & Optimality
Rohan C. Shekhar
Pembroke College
iii
iv
Acknowledements
I would firstly like to thank my supervisor, Professor Jan Maciejowski, for his guidance and
support during my years of research in Cambridge. His knowledge of the field and passion
for his research has been truly inspirational. It has been an honour to work with him as his
student.
Additionally, I would like to thank Dr Edward Hartley for his advice and patience. His
outstanding knowledge of all things MPC has helped me tremendously and I am grateful for
our numerous fruitful discussions, from which emerged many of the ideas presented in this
thesis.
I am also grateful to the Menzies Foundation and the Cambridge Commonwealth Trust for
their financial support throughout my time in Cambridge. I am privileged to be a part of the
community of scholars who have benefitted from these organisations.
Many thanks go to the members of the Control Group for providing an intellectually stimulating
environment for research. I am especially grateful to my fellow students for their friendship
and company throughout the long days prior to submission of this thesis. I wish them well in
their future endeavours.
A special thanks goes to the Pembroke College community. In particular, I would like to
acknowledge all of my friends in the Graduate Parlour, for making my time in Cambridge
truly unforgettable.
Most of all, I would like to thank my family, for their love and unwavering support throughout
the last four years.
v
vi
Declaration
As required by University Statute, I hereby declare that this dissertation is not substantially
the same as any that I have submitted for a degree or diploma or other qualification at any
other university. The dissertation is a result of my own work and includes nothing which is
the outcome of work done in collaboration, except where specified explicitly in the text. I also
declare that the length of this dissertation is less than 65,000 words and the number of figures
is less than 150.
Rohan C. Shekhar
Pembroke College
Cambridge
16th April 2012
vii
viii
Contents
Contents ix
List of Figures xv
List of Tables xvii
Nomenclature xix
1 Introduction 1
1.1 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Notational Conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2 Background 5
2.1 Overview of Model Predictive Control . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.1 Principle of Operation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.2 The Receding Horizon Principle . . . . . . . . . . . . . . . . . . . . . . . 6
2.1.3 Constraint Handling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Mathematical Formulation of MPC . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3 Analysis of Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3.1 Sufficient Conditions from Lyapunov Theory . . . . . . . . . . . . . . . . 10
2.3.2 Linear Quadratic MPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.4 Variable Horizon MPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.4.1 Completion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.4.2 Cost Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4.3 An Alternative Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.5 Robust MPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.5.1 Min-Max Predictive Control . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.5.2 Affine Feedback Policies . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.5.3 Constraint Tightening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.5.4 Tube Model Predictive Control . . . . . . . . . . . . . . . . . . . . . . . . 17
ix
CONTENTS
2.6 Robust VH-MPC with Constraint Tightening . . . . . . . . . . . . . . . . . . . . 18
2.6.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.6.2 Controller Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.6.3 Disturbance Feedback Formulation . . . . . . . . . . . . . . . . . . . . . 21
2.7 Optimisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.7.1 Mixed-Integer Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.7.2 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3 Feasible Contingencies 25
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1.1 Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.2 Feasible Contingencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.2.1 Nominal Control Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.2.2 Robustness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2.3 Linear Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3 Robust MPC with Contingencies . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.3.1 MIP implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.4 Contingency Scenarios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.4.1 Costed Contingencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.4.2 Multiple & Prioritised Contingencies . . . . . . . . . . . . . . . . . . . . 41
3.4.3 State-Dependent Contingencies . . . . . . . . . . . . . . . . . . . . . . . . 42
3.5 Simplifications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.5.1 Fixed Contingency Horizon . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.5.2 Explicit Contingency Feasibility . . . . . . . . . . . . . . . . . . . . . . . 44
3.5.3 Prediction Variable Reduction . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.6 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.6.1 Double-Integrator Vehicle . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.6.2 Cart-Pendulum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.8 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4 Move Blocking for Complexity Reduction 59
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.1.1 Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.2 Problem Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.3 Move Blocking Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
x
CONTENTS
4.4 Blocked VH-MPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.4.1 Nominal Controller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.4.2 Robustness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.5 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.5.1 Input Condensed Prediction Model . . . . . . . . . . . . . . . . . . . . . 74
4.5.2 State & Input Condensed Prediction Model . . . . . . . . . . . . . . . . . 75
4.6 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.6.1 The Partitioned Horizon . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.6.2 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.6.3 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.6.4 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.8 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5 Optimal Constraint Tightening Policies 83
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.1.1 Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.2 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.3 Affine Parameterisation of Tightened Sets . . . . . . . . . . . . . . . . . . . . . . 86
5.4 Maximising the Region of Attraction . . . . . . . . . . . . . . . . . . . . . . . . . 88
5.4.1 Parameterisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
5.4.2 Optimisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
5.5 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
5.5.1 Symmetric Target Region Placement . . . . . . . . . . . . . . . . . . . . . 94
5.5.2 Horizon Length Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . 97
5.5.3 Asymmetric Target Region Placement . . . . . . . . . . . . . . . . . . . . 97
5.5.4 Effect on Cost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
5.5.5 Heuristic Optimisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
5.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
5.6.1 Curse of Dimensionality . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
5.6.2 Computational Complexity . . . . . . . . . . . . . . . . . . . . . . . . . . 102
5.6.3 Volume Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
5.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.8 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
xi
CONTENTS
6 Case Study - Surface Excavation 105
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
6.1.1 Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
6.2 Mechanism Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
6.2.1 Kinematic Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
6.2.2 Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6.2.3 External & Dissipative Forces . . . . . . . . . . . . . . . . . . . . . . . . . 112
6.3 Material Interaction Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
6.3.1 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
6.4 System Identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
6.4.1 Pre-Stabilising Controller . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
6.4.2 Identification of Trajectory Tracking System . . . . . . . . . . . . . . . . . 117
6.5 VH-MPC Controller Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
6.5.1 Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
6.5.2 Cost Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
6.6 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
6.6.1 Normal Operation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
6.6.2 Stall Contingency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
6.6.3 Move Blocking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
6.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
6.8 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
7 Conclusion 129
7.1 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
7.1.1 Feasible Contingencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
7.1.2 Complexity-Reduced VH-MPC with Move Blocking . . . . . . . . . . . 130
7.1.3 Optimal Constraint Tightening Policies . . . . . . . . . . . . . . . . . . . 130
7.1.4 Case Study on Surface Excavation . . . . . . . . . . . . . . . . . . . . . . 130
7.2 Future Research Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
7.2.1 Uncertain Contingencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
7.2.2 Optimised & Multiplexed Move Blocking . . . . . . . . . . . . . . . . . . 131
7.2.3 Cost Optimal Constraint Tightening Policies . . . . . . . . . . . . . . . . 131
7.2.4 Combinations of Strategies . . . . . . . . . . . . . . . . . . . . . . . . . . 131
7.2.5 Shovel & Material Model Generalisation . . . . . . . . . . . . . . . . . . 132
Bibliography 133
xii
CONTENTS
A Constraint Tightening Addenda 139
A.1 Derivation of the Max-Dist Policy . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
A.2 VH-ROA Plots for Different Policies . . . . . . . . . . . . . . . . . . . . . . . . . 140
B Shovel Model & Controller Details 147
B.1 Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
B.2 Controller Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
xiii
xiv
List of Figures
2.1 The receding horizon principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 Qualitative illustration of MPC’s constraint awareness . . . . . . . . . . . . . . . 8
3.1 Illustration of MPC with Contingency Predictions . . . . . . . . . . . . . . . . . 31
3.2 Trajectories with contingency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.3 Primary trajectory velocities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.4 Primary trajectory control inputs . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.5 Primary trajectory fuel use . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.6 Trajectories with mean-costed contingency . . . . . . . . . . . . . . . . . . . . . 50
3.7 Trajectories for two contingency sets . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.8 Trajectories for two contingency sets with priority . . . . . . . . . . . . . . . . . 51
3.9 Trajectories for state-dependent contingency . . . . . . . . . . . . . . . . . . . . 52
3.10 Fuel use for state-dependent contingency . . . . . . . . . . . . . . . . . . . . . . 52
3.11 Diagram of cart-pendulum system . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.12 Step response for cart states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.13 Step response for pendulum states . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.14 Control input to cart . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.1 Values of functions Λ( · , · ) and µ( · , · ) for blocking matrix corresponding to
N = 9 and block lengths {2, 4, 3} . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.2 Illustration of the partitioned horizon regime . . . . . . . . . . . . . . . . . . . . 77
4.3 Cost comparison of different block lengths in the partitioned horizon . . . . . . 79
4.4 Compuation time comparison different block lengths in the partitioned horizon 79
4.5 Closed-loop trajectories for λ = 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.6 Closed-loop trajectories for λ = 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.7 Closed-loop trajectories for λ = 40 . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.8 Control inputs for λ = 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.9 Control inputs for λ = 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.10 Control inputs for λ = 40 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
xv
LIST OF FIGURES
5.1 Four-step ROA for T1 using Pmax-vol . . . . . . . . . . . . . . . . . . . . . . . . . 95
5.2 Total variable horizon ROA for T1 using Pmax-vol . . . . . . . . . . . . . . . . . . 96
5.3 Four-step ROA comparison over different policies for T1 . . . . . . . . . . . . . 96
5.4 Comparison of 4-step ROA for T1 using Pmax-vol with 14-step ROA for T1 using
P2-step-nil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
5.5 Four-step ROA for T2 using Pmax-vol compared with nominal ROA . . . . . . . 98
5.6 Total variable horizon ROA for T2 using Pmax-vol . . . . . . . . . . . . . . . . . . 99
5.7 Four-step ROA comparison over different policies for T2 . . . . . . . . . . . . . 99
6.1 The rope shovel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
6.2 Rope shovel excavation trajectory . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
6.3 Kinematic vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
6.4 Mechanism model showing generalised coordinates . . . . . . . . . . . . . . . . 109
6.5 Static balance of material failure forces . . . . . . . . . . . . . . . . . . . . . . . . 114
6.6 Performance of trajectory tracking controller . . . . . . . . . . . . . . . . . . . . 117
6.7 Performance of identified linear model for crowd length and angle . . . . . . . 118
6.8 Performance of identified linear model for payload . . . . . . . . . . . . . . . . 119
6.9 Actual trajectory and prediction from linear model . . . . . . . . . . . . . . . . . 122
6.10 Crowd position and angle compared to reference . . . . . . . . . . . . . . . . . . 123
6.11 Payload acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
6.12 Hoist and crowd forces from the trim point . . . . . . . . . . . . . . . . . . . . . 124
6.13 Cutting resistance force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
6.14 Trajectory comparison with crowd stall contingency . . . . . . . . . . . . . . . . 125
6.15 Payload comparison with crowd stall contingency . . . . . . . . . . . . . . . . . 125
6.16 Trajectory comparison with move blocking . . . . . . . . . . . . . . . . . . . . . 127
6.17 Payload comparison with move blocking . . . . . . . . . . . . . . . . . . . . . . 127
A.1 Nominal VH-ROA for T1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
A.2 VH-ROA for T1 using max-vol . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
A.3 VH-ROA for T1 using max-dist . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
A.4 VH-ROA for T1 using 3-step nil . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
A.5 VH-ROA for T1 using 2-step nil . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
A.6 Nominal VH-ROA for T2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
A.7 VH-ROA for T2 using max-vol . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
A.8 VH-ROA for T2 using max-dist . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
A.9 VH-ROA for T2 using 3-step-nil . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
A.10 VH-ROA for T2 using 2-step-nil . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
xvi
List of Tables
5.1 ROA volumes over different policies for T1 . . . . . . . . . . . . . . . . . . . . . 95
5.2 ROA volumes over different policies for T2 . . . . . . . . . . . . . . . . . . . . . 98
5.3 Mean costs for Monte Carlo simulation using different policies . . . . . . . . . . 100
5.4 ROA volumes over different policies for 3-state system . . . . . . . . . . . . . . 101
6.1 Kinematic parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
6.2 Dynamic parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
6.3 Material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
xvii
xviii
Nomenclature
Variables, Functions & Constants
0n×m n×m matrix of zeros
1n n-dimensional vector of ones
In n× n identity matrix
δ[i] discrete unit impulse at time i
J∗( · ) value function
Jˆ( · ) suboptimal value function
N prediction horizon length
u input vector
w disturbance vector
x state vector
y output vector
Sets
Bnp(e) n-dimensional p-norm ball with radius e
∅ empty set
QW image of setW under matrix Q ({x | x = Qw, ∀w ∈ W})
N set of natural numbers
N[a,b] set of natural numbers in the range a, . . . , b
R set of real numbers
R+ set of positive real numbers
T terminal / target constraint set
U input constraint set
W state disturbance set
X state constraint set
xix
NOMENCLATURE
Y output constraint set
Z set of integers
Z[a,b] set of integers in the range a, . . . , b
Operators
C matrix horizontal concatenation operator
| · | row-wise absolute value, set cardinality
d · e ceiling function
, defined as
diag{ · } diagonal matrix
diam( · ) polytope diameter
b · c floor function
iff, ⇐⇒ if and only if
=⇒ implies
∩ intersection
⊗ Kronecker product
7→ maps to
‖ · ‖p p-norm
⊕ Minkowski sum (Chapter 3), direct sum (Chapter 4)
	 Pontryagin difference
Pn( · ) projection onto Rn
⊂ subset
⊆ subset or equal to
⊃ superset
⊇ superset or equal to
∪ union
vol( · ) polytope volume
Acronyms
CT Constraint Tightening
DARE Discrete Algebraic Riccati Equation
DLE Discrete Lyapunov Equation
LMI Linear Matrix Inequality
xx
NOMENCLATURE
LP Linear Program
LQR Linear Quadratic Regulator
MILP Mixed-Integer Linear Program
MIP Mixed-Integer Programming
MPC Model Predictive Control
QP Quadratic Program
RMPC Robust Model Predictive Control
ROA Region of Attraction
SDP Semidefinite Program
SOCP Second Order Cone Program
VH-MPC Variable Horizon Model Predictive Control
VH-ROA Variable Horizon Region of Attraction
Definitions
Compact Set: a set which is closed and bounded
Completion: the time at which the system state first enters a specified target set
Move Blocking: an input parameterisation strategy where groups of adjacent-in-time inputs, or
perturbations to a stabilising control law, are constrained to have the same value
Output Admissibility: a property of the state and control input of the system that ensures
satisfaction of output constraints
Polyhedron: a closed, but not necessarily bounded set formed by the intersection of half-spaces
Polytope: a compact set formed by the intersection of half-spaces
Recursive Feasibility: the property of an MPC controller that ensures the existence of a feasible
solution to the optimisation problem at all time steps from an initial feasible solution
Region of Attraction: the set of states for which there exists an initial feasible solution to the
MPC optimisation problem
Subspace: a subset of a vector space that is closed under addition and scalar multiplication
xxi
xxii
Chapter 1
Introduction
Model Predictive Control (MPC) is an optimisation-based control paradigm that has been
used successfully in the process industries for many years at a supervisory level. Its intuitive
formulation and unique constraint-handling ability permit system operation closer to constraint
boundaries, which can improve performance. Recent advances in computational power have
made it possible to implement MPC directly at a low level on systems requiring fast update
rates. This has been further aided by algorithmic advances, which have improved solution
times for the constrained finite-horizon optimisation problem solved at each iteration of the
MPC controller [WB10]. Conditions that guarantee convergence of predictive controllers for
regulation problems are well established in the literature [May+00; RM09]. In such problems
the system’s output is required to asymptotically track some reference value.
However, regulation problems are not always an appropriate description of the desired control
objectives. For applications like vehicle manœuvring, the control objective may be to simply
reach a specified target region, whilst minimising some cost function. Hence, the notion
of stability is replaced with that of completion [RH03]; once the target region is reached (at
completion), the problem ends. In essence, it is the transient that is of importance in problems
of this nature. It is assumed that, if required, another controller is deployed after the completion
time. For most applications of completion problems, it is desirable to have a finite completion
time. Whilst stability is no longer required, it is also desirable to maintain a recursively feasible
optimisation problem at each controller iteration up to completion.
Variable Horizon Model Predictive Control (VH-MPC) [MM93] has been shown to provide
these guarantees. In this form of MPC, the horizon length is a decision variable in the con-
strained optimisation problem solved at each iteration. Recent applications of VH-MPC have
addressed completion problems, where the requirements on terminal set invariance, that are
central to the stability proofs of classical MPC, are relaxed. With Constraint Tightening (CT)
applied for robustness, the guarantees of recursive feasibility and finite-time completion can
be shown to hold for linear systems controlled by VH-MPC even in the presence of bounded
state disturbances [RH03; RH06a].
Interesting research challenges are posed when considering robustness and optimality in real
1
1. INTRODUCTION
world applications of VH-MPC. One problem is ensuring robust finite-time completion under
changing system dynamics and target sets. Such changes to the original control problem could
be caused by faults or changes in operational priorities. When the probable changes are well
defined, it is prudent to ensure the existence of contingency trajectories to the new target sets,
under the new dynamics and constraints. If it is unknown precisely when the control problem
will change, the contingencies need to be robustly available at all times prior to completion of
the original control problem.
Another research opportunity arises when considering the application of robust VH-MPC
to systems requiring prediction over long time horizons. Given that predictions of system
dynamics are required until predicted completion of the control manœuvre, computational
complexity can be an issue for real-time implementation. Techniques like move blocking
[Cag+07] can reduce the computational burden; however, naive application of this technique
will not provide the required robust completion guarantees. New formulations are required to
achieve the aim of reducing complexity and providing the desired theoretical guarantees.
A third research problem presents itself when analysing the effect of robustness modifications
on optimality measures of controller performance. In particular, the manner in which con-
straints are tightened directly affects the size of a VH-MPC controller’s Region of Attraction
(ROA), which is the set of states for which a feasible initial solution exists to the MPC optimisa-
tion problem. In many applications, the size of the ROA is an important measure of control
performance, as it determines the ability of the controller to deal with large ranges of initial
conditions. Tightening the constraints in such a way as to ensure the largest possible ROA size
is desirable.
This thesis develops novel solutions to these three research challenges. It presents a new robust
formulation of VH-MPC with feasible contingencies; develops a robust, flexible move blocking
scheme that preserves recursive feasibility guarantees; and formulates a new algorithm for
constraint tightening that enlarges the controller’s ROA. Numerical examples are used to
illustrate the efficacy of the proposed solutions. A case study is presented to demonstrate
the effect of the techniques on a nontrivial application, namely a model surface excavation
machine.
1.1 Outline
The main body of this dissertation is structured as follows:
Chapter 2 provides a comprehensive overview of the theory that underpins the contributions
in this thesis. It begins with a qualitative overview of model predictive control and intro-
duces classical formulations for regulation problems. The stability properties of these
formulations are also reviewed. Variable horizon completion problems are introduced
next, with an emphasis on how these differ from the classical formulations. An overview
2
1.1. OUTLINE
of popular robustness modifications is then presented, with a detailed outline of con-
straint tightening as applicable to VH-MPC completion problems. The chapter concludes
with details on implementing the variable horizon with mixed-integer programming.
Chapter 3 introduces the notion of feasible contingencies in the context of completion prob-
lems. The sufficient conditions for the design of a controller with nominal and robust
contingency availability are detailed first. Subsequently, a VH-MPC controller is for-
mulated to guarantee contingency availability at all time steps prior to completion.
A modified form of constraint tightening is used to achieve the required robustness
to bounded disturbances during normal operation and under contingency conditions.
The controller formulation is shown to guarantee robust recursive feasibility, robust
finite-time completion and robust contingency availability. Details are also provided for
implementation of the resulting quasi-convex optimisation problems using mixed-integer
programming.
Various contingency scenarios are then presented, including multiple, costed and state-
dependent contingencies. In each of these scenarios, sufficient conditions are derived
for preserving the guarantee of robust finite-time completion. Simplifications to reduce
the computational burden of the optimisation problem are also outlined. The chapter
concludes with numerical examples on two systems that demonstrate the power of the
controller formulation and contingency scenarios.
Chapter 4 presents a technique for reducing the computational complexity of VH-MPC,
whilst still maintaining robust completion guarantees. A new move-blocked formulation
of VH-MPC is developed, which allows the number of input decision variables in the
optimisation problem to be reduced. The variable horizon is achieved by solving multiple
optimisation problems at each time step corresponding to different horizon lengths and
move blocking patterns. A novel constraint tightening methodology is developed to
guarantee robustness to bounded disturbances in the presence of blocking constraints.
Implementational techniques to ensure reductions in computation time with the new
methodology are also outlined.
An example blocking scheme is then proposed, which allows fine-grained control action
during the start and end of a control manœuvre, but applies move blocking in the
middle to reduce complexity. This kind of scheme is suited to vehicle control problems,
in particular, spacecraft rendezvous, where robust constraint satisfaction is potentially
required over long prediction horizons. A simple numerical example is used to illustrate
the effect on computation time and closed loop cost under disturbances of applying the
scheme with varying degrees of complexity reduction.
Chapter 5 presents a way of tightening constraints to improve the region of attraction of a VH-
MPC controller. It shows how the constraint tightening policy, parameterised in terms
of direct disturbance feedback, can be formulated to maximise the volume of an inner
approximation to the controller’s true ROA. Through numerical examples, the technique
is shown to produce vastly improved ROA volumes over other proposed techniques
3
1. INTRODUCTION
for designing constraint tightening policies. It is also shown how the improvement in
ROA volume can drastically decrease the maximum horizon length required to achieve
an initial feasible solution. Simplifications for higher state dimensions and prediction
horizons are also discussed.
Chapter 6 illustrates the application of VH-MPC and the techniques developed in the previous
chapters to a non-trivial system, namely a surface excavation machine. A simplified non-
linear kinematic and dynamic model is developed for a two-degree-of-freedom machine
used widely in surface mining operations. A material model is also developed to simulate
the effect of excavation forces. The machine model is then pre-stabilised with an inner
loop controller, following which a discretised and linearised model of the full machine
and material system is obtained for the application of VH-MPC to trajectory generation.
By appropriately designing the cost function and terminal constraints, realistic digging
trajectories are produced. Further simulations show the effect of including contingency
availability and applying move blocking to the linearised model. The results demonstrate
the applicability of VH-MPC and the theory developed in this thesis to problems other
than vehicle control.
Chapter 7 summarises the conclusions and contributions of the thesis. Future research
directions are also proposed.
Appendix A provides some derivations and additional plots pertaining to Chapter 5.
Appendix B lists the equations of motion and controller parameters for the model develped
in Chapter 6.
1.2 Notational Conventions
Whilst the chapters in this thesis are all related to optimality and robustness in VH-MPC, they
each require different mathematical tools. Hence, adopting a completely coherent system of
notation is impractical. Within each chapter, coherency is preserved; however, some symbols
take on different roles in subsequent chapters. Common nomenclature has been listed at
the start of this thesis. In addition, some of the listed symbols are given a more detailed
introduction at the start of each chapter where appropriate.
4
Chapter 2
Background
This chapter provides an overview of background theory that underlies all of the contributions
in this thesis. To fulfil this aim, it presents a comprehensive review of relevant background
literature. A brief introduction to Model Predictive Control is provided first, following which,
standard formulations of MPC for regulation problems are described and their stability prop-
erties analysed. The notions of a variable horizon and completion are then introduced as an
extension of these formulations.
Robust model predictive control is discussed next, where a brief review is presented of the
most popular techniques for ensuring a priori robustness to bounded disturbances. Emphasis
is given to the constraint tightening approach, which is used extensively in this thesis. In
particular, the formulation of constraint tightening in conjunction with a variable horizon
controller is presented, as applicable to completion problems.
2.1 Overview of Model Predictive Control
This section provides an introduction to Model Predictive Control. A general overview of
its principle of operation is presented first. As part of this overview, the receding horizon
optimisation concept is introduced, which is at the heart of MPC. A qualitative description of
MPC’s constraint handling ability is also given.
2.1.1 Principle of Operation
The basic mechanism of MPC is simple and intuitive. In most common formulations, it is a
discrete-time control methodology, so it only generates control actions at discrete time intervals.
Furthermore, it requires an explicit system model, capable of predicting behaviour over a
number of time steps into the future. This time period is called the prediction horizon. An MPC
controller seeks to find the optimal control input sequence over the prediction horizon that
minimises some desired cost function, using the system model to make predictions about
future behaviour. The salient characteristic of its operation is that it applies the first optimal
5
2. BACKGROUND
input to the system and repeats the process again over the prediction horizon, which is now
shifted forward by a time step. This is the reason that MPC is also known as receding horizon
control. As explained in the survey paper by Mayne et al. [May+00], the essential difference
from standard optimal control is that in general, an optimisation problem is solved at each
time step. This is advantageous, as constraints can be included in the formulation by simply
solving numerically a constrained optimisation problem.
time
in
pu
t Constraint
time
ou
tp
ut Setpoint
Constraint
k k + 1
(a) System at time k
time
in
pu
t Constraint
time
ou
tp
ut Setpoint
Constraint
k k + 1 k + N
Optimal Input
Predicted Output
(b) Prediction made over horizon length N
time
in
pu
t Constraint
time
ou
tp
ut Setpoint
Constraint
k k + 1 k + N
Applied Input
Actual Output
(c) First step of optimal input applied
time
in
pu
t Constraint
time
ou
tp
ut Setpoint
Constraint
k k + 1 k + N
New Optimal Input
New Prediction
(d) Prediction made over receding horizon
Figure 2.1: The receding horizon principle (adapted from [Cam09])
2.1.2 The Receding Horizon Principle
Figure 2.1 illustrates the idea of the receding horizon for a SISO system. The system input
and response up to step k is illustrated in Plot (a). The MPC controller then calculates the
input over the horizon of N steps required to optimise a given cost function, as shown in Plot
(b). Sometimes, the input will be assumed to vary only over the first few steps of the horizon
(known as the control horizon), and then be assumed constant, which can reduce the number
of decision variables over which optimisation is required. Sample operating constraints are
indicated on the figures.
After calculating the optimal input sequence, the predictive controller applies the first element
6
2.2. MATHEMATICAL FORMULATION OF MPC
of this sequence, as shown in Plot (c). In general, the resulting output will be different to the
predicted output, due to disturbances and modelling errors. The notion of a receding horizon
is pictured in Plot (d), where the optimisation process is repeated again to calculate a new
input signal to k + N + 1 steps. The horizon length, being constant, recedes by one time step
after each subsequent control action.
2.1.3 Constraint Handling
Traditional control methods take ad hoc approaches towards handling constraints. Often, this
involves choosing operating points far enough away from constraint boundaries to decrease
the likelihood of violating them. The drawback of such a strategy is that these operating points
may not be optimal. Figure 2.2 illustrates this in a qualitative manner, assuming a linear system
and Gaussian-distributed disturbances.
Graph (a) shows the distribution of outputs for a poorly tuned linear controller. Note that
the mean of the distribution must be kept significantly far from the constraint boundary to
have an acceptably low probability of violation. Graph (b) shows the distribution of outputs
that may be achieved by using an optimal linear controller. The variance is vastly improved,
however the symmetric nature of the distribution still limits the proximity of the operating
point to the constraint boundary. A possible distribution of outputs from an MPC controller is
shown in Graph (c). Note the salient feature of this distribution: it is asymmetric, allowing the
operating point to move closer to the constraint boundary and thereby improving efficiency of
operation, whilst still retaining an acceptably low probability of constraint violation. In effect,
the controller is taking different actions when a disturbance causes the output to approach the
constraint boundary. The resulting control law, however, is nonlinear. Hence, MPC allows
superior handling of constraints at the expense of a nonlinear control law.
2.2 Mathematical Formulation of MPC
This section presents the classical state space formulation of MPC for regulation to the origin,
given a general m-input, n-state and p-output time-invariant system. The derivations that
follow are based on [May+00]. It will be assumed that the system can be approximated by a
system of explicit difference equations
x(k + 1) = f (x(k), u(k)) (2.1)
y(k) = g(x(k)), (2.2)
where f ( · , · ) will be taken to have an equilibrium point at the origin. The vector of states is
denoted x, the input vector is written as u and the function f ( · , · ) defines the successor state
in terms of the current state and input. The function g( · ) represents the system output as a
function of the states alone, assuming there is no input feedforward. Constraints on the input
7
2. BACKGROUND
(a)
(b) (c)
Constraint
Figure 2.2: Qualitative illustration of MPC’s constraint awareness (adapted from [Mac02])
and state vectors are assumed to be of the form
u(k) ∈ U ⊂ Rm (2.3)
x(k) ∈ X ⊂ Rn, (2.4)
where U is usually a convex, compact subset of Rm and X is usually a convex, closed subset of
Rn, where both sets include the origin.
Remark 2.1. In the derivations that follow and indeed in the remainder of this thesis, full state
information is assumed to be available.
The aim of designing a controller is usually to drive the system either to the origin or to an
equilibrium state xr for which g(xr) = r, where r is the desired output value at equilibrium. It
will be assumed without loss of generality that xr = 0, as an affine transformation can be used
to shift the equilibrium state to the origin for a general system. Then, define a prediction model
with dynamics
x(k|k) = x(k) (2.5)
x(k + j + 1|k) = f (x(k + j|k), u(k + j|k)) (2.6)
y(k + j|k) = g(x(k + j|k)), (2.7)
subject to the constraints
u(k + j|k) ∈ U (2.8)
8
2.2. MATHEMATICAL FORMULATION OF MPC
x(k + j|k) ∈ X (2.9)
x(k + N|k) ∈ T , (2.10)
for all j ∈ Z[0,N−1]. The notation x(k + j|k) denotes a prediction j steps into the future of the
state x from the current time k. Constraint (2.10) requires some terminal constraint set T to be
entered at the end of the prediction horizon of length N.
A finite horizon cost function for the prediction model is then typically specified in the form
JN(x(k), u(k)) =
N−1
∑
j=0
`(u(k + j|k), x(k + j|k)) + F(x(k + N|k)), (2.11)
where
u(k) = {u(k|k), u(k + 1|k), . . . , u(k + N − 1|k)} = {u(k + j|k)}N−1j=0 (2.12)
is the sequence of inputs applied over the prediction horizon, x(k) is the current state, `( · , · )
is the stage cost and F( · ) is some terminal cost function.
Remark 2.2. For technical reasons, it is assumed that the stage cost is strictly greater than zero
everywhere except for the origin, where it has the value zero, and decreases with the proximity
of the state to the origin. Mathematically, this requires the existence of a K∞-function c( · ) such
that
`(x, u) ≥ c(‖(x, u)‖). (2.13)
The model predictive control problem can now be formally defined in terms of the cost function.
For the given state of the system at time step k, the aim is to minimise JN( · , · ) whilst satisfying
all constraints on the inputs, intermediate states and terminal state. Given that `( · , · ) and
F( · ) are time invariant, the optimal cost or value function depends only on x(k). This value
function is given by
J∗N(x(k)) = min
u(k)
JN(x(k), u(k))
= min
u(k)
N−1
∑
j=0
`(x(k + j|k), u(k + j|k)) + F(x(k + N|k)), (2.14)
with the minimisation subject to the constraints (2.8) - (2.10). The corresponding input sequence
that minimises the value function is given by the expression
u∗(k) = arg min
u(k)
JN(x(k), u(k))
= {u∗(k + j|k)}N−1j=0 . (2.15)
Note the state prediction at time k + j resulting from applying u∗(k) is denoted x∗(k + j|k).
Only the first optimal input from (2.15) is applied to the system, resulting in the time-invariant
model predictive control law
u(k) = κN(x(k)), (2.16)
9
2. BACKGROUND
where κN( · ) is defined by
κN(x(k)) = u∗(k|k). (2.17)
2.3 Analysis of Stability
Having formally defined the action of predictive control on a general dynamical system, its
stability properties will be analysed. As mentioned in Section 2.1.3, a nonlinear control law
results from MPC with active constraints. This means that finding necessary conditions for
stability is, in general, difficult. Instead, the majority of literature focuses on finding sufficient
conditions for stability through Lyapunov methods. In particular, Mayne et al. [May+00]
outlines the sufficient conditions on the terminal cost F( · ), the terminal constraint set T and a
candidate terminal state feedback control law κ f (x(k)) that guarantee stability. These sufficient
conditions can be derived by utilising the value function J∗N( · ) as a Control Lyapunov Function
(CLF), and placing an upper bound on this function at the following time step (the so-called
“direct” method). In the derivations that follow, the value function will be assumed continuous.
Following the derivations, techniques for satisfying the sufficient conditions for stability will
be discussed.
2.3.1 Sufficient Conditions from Lyapunov Theory
This method of deriving stability conditions seeks to find an upper bound on the value function
at the successor state, namely J∗N(x(k + 1)) = J
∗
N( f (x(k), κN(x(k))). The upper bound is
calculated by defining a feasible control sequence applied at time k + 1 of the form
uˆ(k + 1) = {u∗(k + 1|k), u∗(k + 2|k), . . . , u∗(k + N − 1|k), uˆ(k + N|k + 1)}, (2.18)
which is simply the optimal input sequence (2.15) shifted, with an additional input term as the
new tail. To ensure that the sequence uˆ( · ) is feasible, the tail can be defined by some terminal
control law
uˆ(k + N|k + 1) = κ f (x∗(k + N|k)),
where κ f ( · ) is both input and state admissible. This will be the case provided the terminal
constraint set satisfies
T ⊆ X
and the terminal control law satisfies, for all x ∈ X ,
κ f (x) ∈ U
f
(
x, κ f (x)
) ∈ T .
10
2.3. ANALYSIS OF STABILITY
The resulting state trajectory from applying this sequence is then
xˆ(k + 1) = {x(k + 1), x∗(k + 2|k), . . . , x∗(k + N|k), xˆ(k + N + 1|k + 1)},
where perfect model matching is assumed, so x(k + 1) = x∗(k + 1|k) under the action of MPC.
Then, the new terminal state prediction is given by
xˆ(k + N + 1|k + 1) = f (x∗(k + N|k), κ f (x∗(k + N|k))). (2.19)
The cost at time k + 1 associated with applying the input sequence can be found by subtracting
the old zeroth-stage and terminal costs from J∗N( · ) and then adding the new (N − 1)th stage
and terminal costs, giving
JN(x(k + 1), uˆ(k + 1)) = J∗N(x(k))− `(x(k), κN(x(k)))− F(x∗(k + N|k))
+ `(x∗(k + N|k), κ f (x∗(k + N|k))) + F(xˆ(k + N + 1|k + 1)). (2.20)
As the input sequence uˆ( · ) is feasible, but not necessarily optimal, this cost is an upper bound
for the value function J∗N(x(k + 1)). By requiring the sum of the last three terms of (2.20) to be
less than or equal to zero, closed-loop stability can be ensured. This requires that
F(xˆ(k + N + 1|k + 1))− F(x∗(k + N|k)) + `(x∗(k + N|k), κ f (x∗(k + N|k))) ≤ 0. (2.21)
Observe that (2.21), as well as the technical requirements on `( · , · ) imply that
J∗N(x(k + 1))− J∗N(x(k)) + `(x(k), κN(x(k))) ≤ JN(x(k + 1), uˆ(k + 1))− J∗N(x(k))
+ `(x(k), κN(x(k)))
(2.22)
≤ 0.
Hence, the value function is a CLF if (2.21) is satisfied. All of these requirements are summarised
in Proposition 2.3 [May+00].
Proposition 2.3. If the conditions C1-C4 are satisfied by the terminal set, terminal cost function and
terminal controller, where
C1: T ⊂ X , T is closed and 0 ∈ T (T is state admissible)
C2: κ f (x) ∈ U , ∀x ∈ T (the terminal control law is input admissible for all terminal states)
C3: x+ , f (x, κ f (x)) ∈ T , ∀x ∈ T (T is positively invariant under κ f ( · ))
C4: F(x+)− F(x) + `(x, κ f (x)) ≤ 0, ∀x ∈ T (F( · ) is a local CLF),
asymptotic stability of the origin is ensured.
Remark 2.4. Mathematically, conditions C1-C4 ensure that for all x ∈ X for which an initial
feasible solution to the N-step optimisation problem exists, as k→ ∞, x → 0 in the closed loop.
11
2. BACKGROUND
The fact that given a feasible solution at time k, another one exists at time k + 1 is known as the
recursive feasibility property.
Mayne’s survey paper [May+00] describes the evolution of stabilising modifications to the finite
horizon problem to satisfy C1-C4. The addition of a terminal equality constraint, terminal cost
function, terminal constraint sets or combinations of these have all been shown to theoretically
guarantee stability.
One of the simplest methods is the use of a terminal equality constraint, as outlined by Keerthi
and Gilbert [KG88]. This essentially requires the system to be at the equilibrium state by the end
of the prediction horizon (i.e. x(k + N|k) = 0). In terms of the conditions C1-C4, the terminal
set contains only one element: T = {0}. As 0 ∈ X , C1 is satisfied. The terminal controller is
not required to take any action, since the predicted state trajectory reaches equilibrium in finite
time, so κ f (x(k)) = 0. Hence, C2 is satisfied, given κ f (0) = 0 ∈ U . The third condition C3 is
then also satisfied, given that the equilibrium point at the origin is a positively invariant set
by definition i.e. f (0, 0) = 0 ∈ T . No terminal cost is used in this case, hence F( · ) = 0. In
addition, the technical conditions on `( · , · ) guarantee that `(0, 0) = 0. These ensure that C4
is satisfied: F( f (0, κ f (0)))− F(κ f (0))− `(0, κ f (0)) = 0. It is clear, therefore, that a terminal
equality constraint guarantees stability.
A more general method for guaranteeing stability is known as dual mode MPC [CLM96]. In this
paradigm, the controller imposes a terminal set constraint T satisfying C1 which is invariant
with respect to the terminal control law κ f (x). This ensures satisfaction of C2-C3. The terminal
cost is then chosen to satisfy C4. Once the actual state of the system enters T , the controller is
switched to κ f ( · ), which stabilises the system at the origin. Other exotic variants such as triple
mode MPC have also been proposed [RKC00].
2.3.2 Linear Quadratic MPC
The linear unconstrained case has a particularly nice form. Given the linear system
f (x, u) = Ax + Bu, (2.23)
the stage cost is usually defined by the quadratic form
`(x, u) = xTQx + uTRu, (2.24)
where Q  0 and R  0. Since the system is unconstrained, the conditions C1-C3 are trivially
satisfied. A constant gain terminal controller of the form κ f (x) = K f x is again chosen, where
the gain is selected such that the eigenvalues of A+ BK f lie within the unit circle. This requires
that the pair (A, B) is stabilisable. Then, the terminal cost is chosen to have the form
F(x) = xTPx, (2.25)
12
2.4. VARIABLE HORIZON MPC
where P is the positive definite solution to the Discrete Lyapunov Equation (DLE)
(A + BK f )TP(A + BK f )− P + Q + KTf RK f = 0. (2.26)
Choosing the terminal cost in this way ensures that C4 is satisfied, making the origin asymptot-
ically stable with region of attraction Rn. If the terminal controller is chosen to be
K f = −(BTPB + R)−1BTPA, (2.27)
then (2.26) becomes a Discrete Algebraic Riccati Equation (DARE), with a unique solution for
P, given the conditions on Q and R. The resulting controller then produces the same action as
an infinite horizon LQR controller.
In the linear constrained case, it is explained in [May+00] that for a stable system with no
state constraints and convex, compact input constraints, there is no requirement for a terminal
constraint set, so the reasoning presented above can be used again to prove asymptotic conver-
gence. However, when state and input constraints are enforced, or the system is unstable, a
terminal cost is not sufficient to ensure stability.
To ensure asymptotic convergence under constraints, a dual mode scheme can be used. Define
the terminal constraint set to satisfy C1 and be positively invariant under the linear terminal
control law K f . This requires the terminal set to satisfy
(A + BK f )x ∈ T (2.28)
K f x ∈ U , (2.29)
for all x ∈ T . Then, C1 - C3 are satisfied. By choosing P to satisfy (2.26), C4 is satisfied,
ensuring asymptotic stability of the origin. The terminal controller can be chosen to be the
solution of the infinite horizon unconstrained problem, as in (2.27), which will ensure that the
behaviour of the system once inside the terminal set will be equivalent to the infinite horizon
(unconstrained) optimal solution. With K f chosen in this way, a suitable positively-invariant T
is given by a constraint-admissible level set of xTPx.
2.4 Variable Horizon MPC
Having defined fixed horizon predictive control problems and analysed their stability prop-
erties, the notion of a variable horizon is now introduced. In this variety of MPC, the value
function takes the form
J∗(x(k)) = min
u(k),N(k)
N(k)−1
∑
j=0
`(x(k + j|k), u(k + j|k)), (2.30)
subject to the pointwise-in-time constraints
u(k + j|k) ∈ U (2.31)
13
2. BACKGROUND
x(k + j|k) ∈ X (2.32)
x(k + N(k)|k) ∈ T , (2.33)
for all j ∈ Z[0,N(k)−1], where the horizon length N(k) is a decision variable in the constrained
optimisation problem solved at each k. Essentially, N(k) denotes the first time step at which
the set T is entered. The corresponding optimal horizon length is denoted N∗(k). No terminal
cost is applied in this case.
One of the earliest mentions of a varying prediction horizon is in the paper by Michalska and
Mayne [MM93], where a dual mode scheme is used for control of nonlinear systems. The
observation is made that recursive feasibility is easier to achieve with a variable horizon. In
(2.18), a tail term has to be added to the shifted input sequence to ensure recursive feasibility,
as the horizon length must stay the same. However, allowing the prediction horizon to vary
means that the input sequence
uˆ(k + 1) = {u∗(k + 1|k), u∗(k + 2|k), . . . , u∗(k + N∗(k)− 1|k)} (2.34)
provides a feasible solution to the optimisation problem solved at time k + 1, with a horizon
length N∗(k)− 1.
2.4.1 Completion
Inherent in the definition of the variable horizon is the idea of completion, which denotes
that the state of the system is driven into the terminal set T in N(k) steps. This can form the
first part of a dual mode scheme, where the control action will be switched once the target set
is entered to stabilise at an equilibrium point. It can be shown that with an initially feasible
solution and the stage cost defined to satisfy (2.13), completion is achieved in finite time. The
technical requirement (2.13) ensures that for all j < N∗(k), there exists some constant c ∈ R+
such that
`(x(k + j|k), u(k + j|k)) ≥ c, (2.35)
as the state prediction x(k + j|k) is outside the terminal set. It follows that, using the shifted
feasible input at time k + 1 with horizon length N∗(k)− 1,
J∗(x(k + 1))− J∗(x(k)) ≤ Jˆ(x(k + 1))− J∗(x(k))
= −`(x(k), u∗(k|k))
≤ −c. (2.36)
This means that the value function at each time step decreases by at least c > 0. As the value
function is non-negative by construction, this guarantees that the set T must be entered in at
most bJ∗(x(k))/cc steps, from current state x(k). Note that this guarantee does not require the
satisfaction of C1-C4, which pertain to stability.
Remark 2.5. The value function in this case is discontinuous, as for any state within the
14
2.4. VARIABLE HORIZON MPC
terminal state, the horizon length, and therefore the value function, is zero. For this reason,
J∗( · ) will be referred to as a Lyapunov-like function when the variable horizon is used.
2.4.2 Cost Functions
For a stage cost chosen to satisfy (2.13), if an initial feasible solution exists, completion is
achieved in finite time. However, this stage cost has been designed with the objective of
stabilising at the origin. In some applications, there is no requirement for stability; enforcing
the invariance requirement on T (C3) can be overly restrictive for such applications [RH03;
RH06a]. In vehicle manœuvring problems, the control objective is often to simply achieve
completion in some target set T , with the control problem ending after this occurs. This thesis
is concerned mainly with control problems of this form.
The previous subsection showed that the key requirement for finite-time completion is a
positive lower bound on the stage cost outside T . Hence, different kinds of cost functions
can be defined, which still guarantee finite-time completion, but have no requirement on T
containing the origin. Since the horizon length is a variable in the optimisation problem, a
simple minimum-time value function [SM98] of the form
J∗(x(k)) = min
N(k)
N(k), (2.37)
can be defined, which will optimise the time to completion. Of course, this may result in
aggressive control action, so the cost function can be modified to penalise a weighted sum of
the control inputs over the prediction and the time to completion. The value function then
takes the form
J∗(x(k)) = min
u(k),N(k)
N(k)−1
∑
j=0
1+ γ ‖u(k + j|k)‖p , (2.38)
for some p-norm of the inputs weighted by γ > 0. In [RH03; RH06a], a 1-norm is chosen, as
the application in question is spacecraft rendezvous, where this quantity provides a measure
of fuel use during a control manœuvre. Note that no terminal cost has been included, as the
control problem ends at completion.
Remark 2.6. From here on, the p subscript on the norm will be dropped unless a particular
norm is implied.
It is easy to see that when the time to entry into T is included in the cost function, finite-time
completion is guaranteed. From the feasible input sequence at time k + 1 (2.34), the cost
function is guaranteed to reduce by at least unity, which guarantees finite-time entry into T ,
given that the cost function is non-zero by construction. Another important observation is
that these new cost functions make it possible to define the completion problem directly in
terms of the original target set ,without requiring a change of coordinates to have the origin
contained within T . The benefits of being able to do this will become apparent in Chapter 3,
where multiple target sets are defined within the same MPC problem.
15
2. BACKGROUND
Remark 2.7. Note that the new cost functions do not prevent VH-MPC being used for regula-
tion problems. By appropriately choosing T as a control invariant set, a dual mode scheme can
still be used to guarantee asymptotic convergence to a given point within the terminal set.
2.4.3 An Alternative Formulation
The dual mode schemes detailed in [MM93] and [SMR99] also use a variable horizon, albeit
formulated in a slightly different way. A fixed horizon length is used for the optimisation;
however, the schemes multiply the original stage cost by the indicator function α( · ) for the
state being outside the terminal set, giving the new stage cost
¯`(x, u) = α(x)`(x, u), (2.39)
where
α(x) =
0 x ∈ T1 otherwise. (2.40)
This ensures that for a regulation problem, C4 is satisfied and the system is stabilised. Es-
sentially, the trajectory is only costed whilst the state is outside the terminal set, so this is
equivalent to applying costs over a variable horizon up to the completion time. The main
difference is that constraints are still applied after this time, so if a completion problem is
considered with no subsequent stabilisation, the closed loop behaviour may differ from a
controller that explicitly optimises over the horizon length.
2.5 Robust MPC
The preceding sections have presented a formulation of MPC where no a priori anticipation of
disturbances is incorporated. Hence, simply using a shifted, truncated input sequence will not
guarantee recursive feasibility under the action of disturbances. This is due to the fact that at
the next prediction step, the system state will not necessarily be at its one-step prediction from
the current time. To address this, many forms of Robust Model Predictive Control (RMPC)
have been proposed. This section summarises the main contributions, for an additive state
disturbance w(k) ∈ W ⊂ Rn, where the disturbance setW is compact and contains the origin.
2.5.1 Min-Max Predictive Control
Min-max robustness modifications centre around the idea of minimising the cost that would
be incurred by the worst case disturbance realisation over the prediction horizon. By finding
an input sequence that achieves this objective and results in a constraint-admissible trajectory
for all disturbance realisations, robust feasibility is ensured by construction. This so-called
open-loop approach can be excessively conservative, hence a modification that optimises over
16
2.5. ROBUST MPC
state feedback policies instead is outlined in [SM98], and subsequently extended in [KM04].
Essentially, predicted trajectory sequences are formulated for every disturbance realisation,
which is made tractable ifW is a polytope. The variable horizon case is considered at the end
of [SM98]. An alternative approach using Linear Matrix Inequalities (LMIs) is described by
[KBM96], where the pair (A, B) is taken to lie within a known polytope or an ellipsoid. An
upper bound on the cost function is then minimised over state feedback policies.
2.5.2 Affine Feedback Policies
Optimisation over affine feedback policies is described in [GKM06]. This approach differs
from the min-max strategies, as it includes an extra affine term in the control policy at each
prediction step. The set of affine state feedback policies is non-convex; however, it is shown
that the set of affine disturbance feedback policies is equivalent, which results in a convex
optimisation problem. These policies are shown to guarantee robust constraint satisfaction for
bounded disturbances when applied in a receding horizon manner.
2.5.3 Constraint Tightening
This robustness modification for linear time-invariant systems was first elucidated by [GKR97]
and subsequently developed by [CRZ01] and [RH03; Ric05b; RH06a]. A formulation for LTV
systems is described in [Ric05a], and a disturbance feedback parameterisation is introduced by
[KRH07]. In all cases, the essence of Constraint Tightening (CT) is that constraint sets at future
steps into the prediction horizon are shrunk, to allow enough control authority to correct for
the action of an unknown but bounded state disturbance. This can be seen as a partitioning of
the control authority; one part is used to steer the system state into T and the remainder is held
in reserve to mitigate the effect of the error between the predicted and actual state evolution
at each step. By designing T as a robust control invariant set, robust asymptotic convergence
can be guaranteed [RH06b]. For completion problems, robust finite-time completion can be
guaranteed for an arbitrary choice of closed terminal set [RH06a]. A detailed overview of
constraint tightening for linear systems will be presented in the following section.
Remark 2.8. The disturbance feedback formulation of CT is identical to applying the affine
feedback formulation of [GKM06] with time invariant disturbance feedback gains.
2.5.4 Tube Model Predictive Control
This form of robustness modification is closely related to constraint tightening. The basic
formulation, known as disturbance invariant tube MPC, is outlined in [ML01]. It involves
partitioning each term of the input sequence into a feedforward term and a feedback term.
Specifically, the input to the system at time k takes the form
u(k) = u∗(k|k) + K(x(k)− x∗(k|k)), (2.41)
17
2. BACKGROUND
where the constraint x∗(k|k) = x(k) specified in (2.5) is relaxed, and K is a stabilising controller
for (A, B). The relaxation means that the predicted states trajectory now defines a “tube centre”,
which does not necessarily need to originate at the current state. This nominal predicted
trajectory is constrained to reach the terminal set in N steps under the predicted input sequence
u∗(k). The additional feedback term regulates the true state of the system to the predicted tube
centre.
To ensure that u(k) is feasible, the constraints applied to the nominal trajectory need to be
tightened. Generalisations to the basic formulation can be found in [Lan+04], where the
variable horizon case is addressed. Tube MPC has also been applied to nonlinear systems, as
described in [May+11]. In the case of linear systems, CT with a fixed state feedback gain is in
general, less conservative than tube MPC [Tro09].
Remark 2.9. Observe that the former two robustness modifications presented in this section
both increase online computational complexity, whereas the latter two do not.
2.6 Robust VH-MPC with Constraint Tightening
In this section, an n-state, m-input linear system of the form
x(k + 1) = Ax(k) + Bu(k) + w(k) (2.42)
will be assumed, where
w(k) ∈ W (2.43)
is a state disturbance, known to lie within the compact set W , which contains the origin.
Constraints are now generalised to be in terms of an output vector
y(k) = Cx(k) + Du(k), (2.44)
where C and D are a design choice. Hence, the system is subject to constraints of the form
y(k) ∈ Y , (2.45)
where Y is closed. In the remainder of this thesis, whenever a linear system is introduced with
outputs y( · ), it is assumed that these have been designed for the application of constraints.
They do not represent measurements, as full state information is assumed.
Remark 2.10. Individual state and input constraints can be easily cast into output constraint
form. To see this, given the separate state and input constraint sets defined in (2.3) - (2.4), it is
clear that
(x, u) ∈ X × U . (2.46)
18
2.6. ROBUST VH-MPC WITH CONSTRAINT TIGHTENING
By choosing
C =
[
In 0m×n
]T
(2.47)
D =
[
0n×m Im
]T
(2.48)
Y = X × U , (2.49)
the constraint (2.46) can be written in the form (2.45). The output constrained form is more
general, as it permits cross-constraints between states and inputs, which could not be achieved
with individual state and input constraint sets.
2.6.1 Formulation
Using a similar notation to [RH06a], define the matrices
L(0) = In (2.50a)
L(j + 1) = (A + BK(j))L(j) (2.50b)
Q(j) = (C + DK(j))L(j) (2.50c)
and the corresponding sets
Y(0) = Y (2.51a)
T (0) = T (2.51b)
Y(j + 1) = Y(j)	Q(j)W (2.51c)
T (j + 1) = T (j)	 L(j)W . (2.51d)
The sequence of matrices
{K(j) ∈ Rn×m}∞j=0, (2.52)
is known as the state feedback constraint tightening policy. In [RH06a], nilpotent policies are
suggested, as they minimise the number of tightened sets to be calculated. For a degree-h
nilpotent policy, L(h) = 0, so there is no additional tightening applied for j > h.
The prediction model is then defined by the dynamics
x(k + j + 1|k) = Ax(k + j|k) + Bu(k + j|k) (2.53a)
y(k + j|k) = Cx(k + j|k) + Du(k + j|k), (2.53b)
subject to the tightened constraints
y(k + j|k) ∈ Y(j) (2.54a)
x(k + N(k)|k) ∈ T (N(k)). (2.54b)
19
2. BACKGROUND
The operator ‘	’ denotes the Pontryagin set difference, or P-difference, which is also known as
the set erosion operation used in mathematical morphology [Rak+06]. It is defined for sets A
and B by
A	B = {x | x + b ∈ A, ∀b ∈ B}. (2.55)
The properties of the P-difference are described in [KG95]. Implementation of the operator for
polyhedral sets is discussed in [Ker00]. The robust model predictive controller now minimises
a cost function of the form (2.38), subject to the dynamics (2.53) and the operating constraints
(2.54).
2.6.2 Controller Properties
The following theorems describe the recursive feasibility and finite-time completion guarantees
provided by this constraint tightened controller.
Theorem 2.11 (Robust Recursive Feasibility). Given the optimal input, state and output sequences
to the system at time k, namely {u∗(k + j|k)}N∗(k)−1j=0 , {x∗(k + j|k)}N
∗(k)
j=0 and {y∗(k + j|k)}N
∗(k)−1
j=0 ,
with optimal horizon length N∗(k), a candidate suboptimal, but feasible trajectory at time k + 1 is
defined by the relationships
uˆ(k + j + 1|k + 1) = u∗(k + j + 1|k) + K(j)L(j)w(k), ∀j ∈ Z[0,N∗(k)−2] (2.56a)
xˆ(k + j + 1|k + 1) = x∗(k + j + 1|k) + L(j)w(k), ∀j ∈ Z[0,N∗(k)−1] (2.56b)
yˆ(k + j + 1|k + 1) = u∗(k + j + 1|k) + (C + DK(j))L(j)w(k), ∀j ∈ Z[0,N∗(k)−2] (2.56c)
Nˆ(k + 1) = N∗(k)− 1. (2.56d)
Proof. The full proof is provided in [RH06a]. It shows that the candidate solution satisfies
the dynamics constraints, and is admissible with respect to the tightened constraints. Then,
by induction, this implies that the optimisation is feasible until completion, given an initial
feasible solution. 
Theorem 2.12 (Robust Finite-Time Completion). If a feasible solution exists at time k and the
weighting on the cost function (2.38) satisfies
0 < η = 1− γ
(
max
w∈W
∞
∑
j=0
‖K(j)L(j)w‖
)
, (2.57)
then the target set T will be reached in at most bJ∗(x(k))/ηc steps, from the current state.
Proof. See [RH06a]. The cost of the candidate solution (2.56) is compared to the optimal cost at
time k. Using the triangle inequality, this places a lower bound on the amount by which the
cost decreases as in (2.36). Treating the cost as a Lyapunov-like function, if the property (2.57)
is satisfied, the set T must be entered in finite time as J∗( · ) is non-negative by construction
and takes on the value zero within T . The maximum number of steps to completion is then
given by bJ∗(x(k))/ηc. 
20
2.7. OPTIMISATION
Remark 2.13. A p-norm cost is used rather than a quadratic cost, as the lower bound on the
value function under the action of disturbances in Theorem 2.12 requires the stage cost to
satisfy the triangle inequality. The quadratic cost, being a squared norm rather than a norm,
does not satisfy this property.
Corollary 2.14 (Lower Bound on Value Function). For any state x /∈ T for which an initial feasible
solution exists,
J∗(x) ≥ η (2.58)
Proof. The result follows from the fact that the value function must reduce by at least η at each
step and is zero within the terminal set. 
2.6.3 Disturbance Feedback Formulation
In [KRH07], a constraint tightening policy parameterised in terms of direct feedback on the
disturbance is presented. It is shown that all state feedback tightening policies can be represen-
ted in this form, but not all disturbance feedback policies have an equivalent state feedback
representation. Furthermore, it is shown that the expressions for the tightening matrices are
now convex functions of the disturbance feedback policy. For these reasons, the disturbance
feedback approach is used exclusively in this thesis. The detailed formulation is introduced in
subsequent chapters, so is omitted here.
2.7 Optimisation
It can be shown that the standard linear-quadratic fixed horizon MPC problem can be expressed
in the form of a Quadratic Program (QP), if the constraint sets are formed by the intersection
of linear half-spaces [Mac02]. Similarly, a 1-norm cost function will result in a linear program.
However, the variable horizon introduces an additional complication, as the cost function is
not, in general, a convex function of the horizon length. Two main methods are employed to
handle this non-convexity. It is possible to simply evaluate the fixed horizon value function
corresponding to each horizon length and choose the one with a minimum cost. Alternatively,
Mixed-Integer Programming (MIP) may be used to solve a quasi-convex optimisation problem
[RH06a; RH03; BM99]. This section provides an overview of the latter. To ensure a tractable
formulation, a maximum horizon length needs to be specified, which will be denoted N¯.
2.7.1 Mixed-Integer Formulation
The prediction model is augmented with the binary variable sequence
b(k) = {b(k + j|k) ∈ {0, 1}}N¯j=0, (2.59)
21
2. BACKGROUND
where the variable b(k + j|k) is set to unity when completion is predicted at step j. Then, the
VH-MPC optimisation problem can be defined as solving
J∗(x(k)) = min
u(k),b(k)
N¯−1
∑
j=0
jb(k + j|k) + γ ‖u(k + j|k)‖ , (2.60)
subject to the tightened mixed logic and dynamics constraints
j
∑
l=0
b(k + l|k) = 0 =⇒ y(k + j|k) ∈ Y(j) (2.61a)
b(k + j|k) = 1 =⇒ x(k + j|k) ∈ T (j) (2.61b)
N¯
∑
l=0
b(k + l|k) = 1. (2.61c)
Constraint (2.61a) enforces the output constraints only whilst completion has not been reached,
constraint (2.61b) ensures that the binary variable is set to unity only when the state is within
the terminal set and (2.61c) mandates that completion is achieved within N¯ steps.
Remark 2.15. The constraint (2.61c) imposes a particular form on the optimal solution, namely
that completion is reached at precisely one step along the prediction horizon. However, this
does not require the state to leave the terminal set after this time, as (2.61b) is only a one-way
implication. Furthermore, if the (2.61c) is relaxed to allow the sum of the binary variables to be
greater than or equal to unity, the trajectory that results from solving (2.60) will still have the
same form by optimality.
2.7.2 Implementation
If the output and terminal constraints are in half-space form, a “big-M” formulation can be used
to implement the logical implications [Sch+01]. If the output constraint sets can be expressed
in the form
Y(j) = {y | E(j)y ≤ f (j)}, (2.62)
for matrices E( · ) and vectors f ( · ), then (2.61a) can be written as
E(j)y(k + j|k) ≤ f (j) + M
(
j
∑
l=0
b(k + l|k)
)
1, (2.63)
for some suitably large number M. This ensures relaxation of the constraint once completion is
reached. Similarly, if the terminal constraint sets can be expressed in the form
T (j) = {x | G(j)x ≤ h(j)}, (2.64)
then (2.61b) can be written in the form
G(j)x(k + j|k) ≤ h(j) + M(1− b(k + j|k))1, (2.65)
22
2.7. OPTIMISATION
which ensures that if b(k + j|k) is set to unity, the appropriate terminal constraint must be
satisfied.
With a 1-norm objective function (2.38), the resulting optimisation problem solved at each
iteration is a Mixed-Integer Linear Program (MILP). If a 2-norm is used instead, the resulting
problem is a Mixed-Integer Second Order Cone Program (MISOCP). These optimisation meth-
ods rely on branch-and-bound techniques [Sch86] to intelligently solve the convex problems
corresponding to fixed sequences of values for b(k).
23
24
Chapter 3
Feasible Contingencies
3.1 Introduction
In real-world systems, it is important to consider the robustness of a controller to changes in
the control objective, system dynamics, constraints, or a combination of these. In all cases, it is
desirable to have one or more contingencies available, if the changing objectives, dynamics and
constraints are well defined. An example of such a scenario is fault tolerance, where system
parameters under fault may be known, but the exact time at which the fault will occur is not.
If and when the changes occur, the appropriate contingencies can then be activated, that is, a
new control problem can be posed under the changed conditions. The new problem will only
be feasible, however, if the system state at the time of contingency activation lies within the
modified region of attraction.
It is therefore prudent to consider the design of a controller that ensures the availablilty of
feasible contingencies at all times prior to completion of the original control problem. For
practical applicability, it is also important for the controller to be robust to disturbances.
To design such a controller, this chapter builds upon ideas from prior work on nominal
MPC formulations having additional safety constraints. Passive and active fault-tolerance
formulations for spacecraft rendezvous that guarantee collision free abort trajectories in case
of thruster failure are presented in [BH08] and [BH07] respectively. The strategy detailed in
[SHF04], on the other hand, ensures the existence of a safe periodic trajectory at the end of
every prediction horizon, as applied to vehicle path planning. A third idea is presented in
[Car+08] for a continuous time predictive controller with a reactive safety mode, where the
ability to track some safety reference trajectory is always available. Robustness is provided
with an adapted tube formulation.
In combining elements of these strategies, the formulation outlined in the following sections
extends the idea of fault-tolerant escape trajectories to more general contingency scenarios,
by ensuring that the system state can be steered to one or more well defined target regions
if contingency activation is required. It does so by explicitly guaranteeing the existence of
controlled trajectories, or feasible contingencies, at every future state prediction. The approach
25
3. FEASIBLE CONTINGENCIES
in this chapter differs from the active and passive abort strategies by using receding horizon op-
timisation as opposed to single-shot optimisation for path planning. This allows robustness to
bounded disturbances to be provided through a modified constraint tightening approach. Vari-
ous scenarios, including costed, prioritised and state-dependent contingencies are considered.
Two example systems are then used to demonstrate the efficacy of the formulation.
3.1.1 Nomenclature
In this chapter, the symbol ‘⊕’ denotes the Minkowski sum, also known as the dilation
operation in mathematical morphology. It is defined for sets A and B by
A⊕B = {x | x = a + b, ∀a ∈ A, b ∈ B}. (3.1)
3.2 Feasible Contingencies
This section poses the control problem for a general discrete-time system and states the condi-
tions on a control law under which a contingency is guaranteed to be available, for the nominal
and robust cases respectively.
3.2.1 Nominal Control Problem
Consider a switched discrete-time system of the form
x(k + 1) =
 f (x(k), u(k)) k < ksf˜ (α)(x(k), u(k)) k ≥ ks, (3.2)
y(k) =
g(x(k), u(k)) k < ksg˜(α)(x(k), u(k)) k ≥ ks, (3.3)
having state vector x( · ) ∈ Rn, input vector u( · ) ∈ Rm and output vector y( · ) ∈ Rp, where ks
is some unknown switching time and α ∈ N[1,nc] is an unknown parameter, assuming there
are nc possible contingencies that may need to be activated. The output vector is subject to the
pointwise-in-time switched constraints
y(k) ∈
Y k < ksY˜ (α) k ≥ ks, (3.4)
where Y and Y˜ (α) are closed, for all α. The primary control problem is to find an output-
admissible feedback policy u(k) = κ(x(k)) that drives the system state to some target set T in
finite time kc ≤ N¯, given a maximum horizon length of N¯ steps. If the dynamics switch and the
contingency is activated before T has been reached, that is, if ks < kc, the new control problem
26
3.2. FEASIBLE CONTINGENCIES
is to drive the state to a contingency target set T˜ (α) in at most M¯(α) steps after contingency
activation, with the new dynamics, constraints and disturbance input specified by the particular
value of α.
Remark 3.1. A technical assumption is imposed on the sets Y and Y˜ (α), namely that the
intersection of the output admissible sets of states corresponding to both sets is non-empty.
This ensures that there exists a feasible region within which the state at the time of contingency
activation is compatible with both sets of constraints.
Remark 3.2. In the control problem, the dynamics and constraints are assumed to only switch
once, at the time ks.
For notational convenience, define the input vector µ(i), state vector φ(i) and output vector
ψ(i), for i ≥ 0 as
µ(i) = u(ks + i) (3.5)
φ(i) = x(ks + i) (3.6)
ψ(i) = y(ks + i), (3.7)
which will be used to represent the dynamics of the system after the activation of contingency
α. The non-negativity of i implies that these variables only have meaning after contingency
activation. Using this notation, define the q-step controllable set [Ker00]
Kq(T˜ (α)) =
{
φ(0) ∈ Rn ∣∣ ∃{µ(i)}q−10 : {ψ(i) ∈ Y˜ (α)}q−10 , φ(q) ∈ T˜ (α)}. (3.8)
This is the set of all states for which there exists an admissible input sequence after activation of
contingency α that guides the system state admissibly to some target set T˜ (α) in exactly q steps.
Definition 3.3 (Nominal Contingency Feasibility). For the contingency α to be nominally
feasible, the set
FM¯(T˜ (α)) =
M¯(α)⋃
q=1
Kq(T˜ (α)) (3.9)
must be nonempty.
Criterion 3.4 (Nominal Contingency Availability). For a contingency to be nominally available,
it must be nominally feasible and the control input u(k) applied to the system must ensure that
x(k) ∈ FM¯(α)(T˜ (α)) (3.10)
for all 0 < k < kc.
Remark 3.5. If there is no time limit on contingency completion, it is possible to relax the
completion time constraint by letting M¯(α) → ∞ in (3.9), if the limit exists.
27
3. FEASIBLE CONTINGENCIES
3.2.2 Robustness
A bounded state disturbance is now introduced into the formulation. Using a slight abuse of
notation, the state dynamics are altered to
x(k + 1) =
 f (x(k), u(k), w(k)) k < ksf˜ (α)(x(k), u(k), w(k)) k ≥ ks, (3.11)
where the disturbance satisfies the pointwise-in-time switched constraints
w(k) ∈
W k < ksW˜ (α) k ≥ ks. (3.12)
The control problem is now to find an admissible feedback policy u(k) = κ(x(k)) that drives
the system state to the target set T in finite time kc, whilst allowing the state to be driven to the
set T˜ (α) in M¯(α) steps or fewer, for any allowable disturbance realisation. Again, for notational
convenience, define the disturbance after the dynamics switch as
ω(i) = w(ks + i). (3.13)
Now define the robust q-step controllable set [Ker00]
K˜q(T˜ (α)) =
{
φ(0) ∈ Rn ∣∣ ∃{µ(i) = κ˜(φ(j))}q−10 :
{ψ(i) ∈ Y˜ (α)}q−10 , φ(q) ∈ T˜ (α), ∀ω(i) ∈ W˜ (α)
}
. (3.14)
This is the set of all states for which there exists a control law κ˜( · ) after activation of contingency
α that guides the system state admissibly to the target set T˜ (α) in exactly q steps.
Definition 3.6 (Robust Contingency Feasibility). For the contingency α to be robustly feasible,
the set
F˜M¯(α)(T˜ (α)) =
M¯(α)⋃
q=1
K˜q(T˜ (α)) (3.15)
must be nonempty.
Criterion 3.7 (Robust Contingency Availability). For the contingency α to be robustly available,
it must be robustly feasible and the control law applied in the closed loop u(k) = κ(x(k)) must
ensure that for all 0 < k < kc,
x(k) ∈ F˜M¯(α)(T˜ (α)). (3.16)
3.2.3 Linear Case
Having introduced the feasible contingency criteria for a general system, the particular case of
a linear system will be analysed. The dynamics prior to contingency activation are described
28
3.3. ROBUST MPC WITH CONTINGENCIES
by the equations
f (x(k), u(k)) = Ax(k) + Bu(k) + w(k) (3.17)
g(x(k), u(k)) = Cx(k) + Du(k), (3.18)
where A ∈ Rn×n, B ∈ Rn×m, C ∈ Rp×n, and D ∈ Rp×m. The ideal control objective is to
admissibly reach the set T in finite time, from the initial state x(0), whilst minimising a p-norm
cost function of the form
J(x(0), y) =
kc−1
∑
k=0
(1+ ‖Θ(y(k)− yr)‖), (3.19)
where y = {y(0), y(1), . . . , y(N − 1)}, Θ ∈ Rp×p and yr is some reference output, which is not
required to correspond to an equilibrium state. It is required that
x(kc) ∈ T . (3.20)
The cost function J( · , · ) penalises a weighted sum of output-dependent operating costs
and time to completion. Using this form of cost function allows more general costs to be
encapsulated than 2.38.
Assuming linear time-invariant dynamics, the system after switching can be similarly repres-
ented by the equations
f˜ (α)(x(k), u(k)) = A˜(α)x(k) + B˜(α)u(k) + w(k)
g˜(α)(x(k), u(k)) = C˜(α)x(k) + D˜(α)u(k),
(3.21)
where A˜(α) ∈ Rn×n, B˜(α) ∈ Rn×m, C˜(α) ∈ Rp×n, and D˜(α) ∈ Rp×m.
Remark 3.8. It is assumed that the pairs (A, B) and (A˜(α), B˜(α)) are controllable.
3.3 Robust MPC with Contingencies
An MPC controller will now be formulated to ensure reachability of the target set T˜ (α) within
M¯(α) steps to satisfy Criterion 3.7. Only the robust case is considered as the nominal problem
is a special case with w(k) = 0 for all k. Define the nominal prediction model under normal
operation by
x(k|k) = x(k) (3.22)
x(k + j + 1|k) = Ax(k + j|k) + Bu(k + j|k) (3.23)
y(k + j|k) = Cx(k + j|k) + Du(k + j|k). (3.24)
29
3. FEASIBLE CONTINGENCIES
It is assumed that the maximum prediction horizon is N¯. Also define the prediction model
under a change in dynamics at prediction step j ≥ 1 by the equations
φ
(α)
j (0|k) = x(k + j|k) (3.25)
φ
(α)
j (i + 1|k) = A˜(α)φ(α)j (i|k) + B˜(α)µ(α)j (i|k) (3.26)
ψ
(α)
j (i|k) = C˜(α)φ(α)j (i|k) + D˜(α)µ(α)j (i|k). (3.27)
These variables predict the behaviour of the system for one particular contingency α.
Recursively define the matrices
L(0) = In (3.28)
L(j + 1) = AL(j) + BP(j) (3.29)
Q(j) = CL(j) + DP(j) (3.30)
for the primary trajectory and
L˜(α)(0) = In (3.31)
L˜(α)(i + 1) = A˜(α) L˜(α)(i) + B˜(α)P˜(α)(i) (3.32)
Q˜(α)(i) = C˜(α) L˜(α)(i) + D˜(α)P˜(α)(i) (3.33)
for contingency α, as well as the corresponding sets
Y(0) = Y , Y(j + 1) = Y(j)	Q(j)W (3.34)
T (0) = T , T (j + 1) = T (j)	 L(j)W (3.35)
for the primary trajectory and
Y˜ (α)0 (0) = Y˜ (α), Y˜ (α)0 (i + 1) = Y˜ (α)0 (i)	 Q˜(α)(i)W˜ (α) (3.36)
T˜ (α)0 (0) = T˜ (α), T˜ (α)0 (i + 1) = T˜ (α)0 (i)	 L˜(α)(i)W˜ (α) (3.37)
Y˜ (α)j+1(i) = Y˜ (α)j (i)	 Q˜(α)(i)L(j)W (3.38)
T˜ (α)j+1 (i) = T˜ (α)j (i)	 L˜(α)(i)L(j)W (3.39)
for the contingency, where j ∈ Z[0,N¯] and i ∈ Z[0,M¯(α)].
In a similar manner to the state-feedback policy introduced in Section 2.6, the sequence of
matrices {P(j)} defines the constraint-tightening policy in disturbance-feedback form. Addi-
tionally, the sequence {P˜(α)(i)} denotes the policy corresponding to contingency α. Tightening
constraints to ensure that these policies can be feasibly implemented guarantees robustness to
the bounded disturbances before and after contingency activation.
30
3.3. ROBUST MPC WITH CONTINGENCIES
The MPC optimisation problem is then to solve
J∗(x(k)) = min
θ
N(k)−1
∑
j=0
(1+ ‖Θ(y(k + j|k)− yr)‖) (3.40)
subject to
y(k + j|k) ∈ Y(j) (3.41a)
x(k + N(k)|k) ∈ T (N(k)) (3.41b)
0 ≤ N(k) ≤ N¯ (3.41c)
ψ
(α)
j (i|k) ∈ Y˜ (α)j (i) (3.41d)
φ
(α)
j (M
(α)
j (k)|k) ∈ T˜ (α)j (M(α)j (k)) (3.41e)
0 ≤ M(α)j (k) ≤ M¯(α), (3.41f)
where j ∈ Z[0,N(k)−1] for the primary prediction variables, j ∈N[1,N(k)−1] for the contingency
prediction variables, i ∈ Z
[0,M(α)j (k)−1]
and Θ ∈ Rp×p is a weighting matrix. The optimisation is
with respect to the set of decision variables,
θ =
⋃
i,j
{
u(k + j|k), µj(i|k), N(k), M(α)j (k)
}
, (3.42)
where it is assumed that nc = 1. Figure 3.1 illustrates the relationships between the state
prediction variables, where the solid arrows indicate dynamics under normal operation and
the dashed arrows indicate the dynamics of contingency α.
x(k)
x(k + 1|k) φ(α)1 (1|k) φ(α)1 (2|k) · · · φ
(α)
1
(
M(α)1 (k)
∣∣∣ k)
x(k + 2|k) φ(α)2 (1|k) φ(α)2 (2|k) · · · φ
(α)
2
(
M(α)2 (k)
∣∣∣ k)
...
x(k + N(k)− 1|k) φ(α)N(k)−1(1|k) φ
(α)
N(k)−1(2|k) · · · φ
(α)
N(k)−1
(
M(α)N(k)−1(k)
∣∣∣ k)
x(k + N(k)|k)
Figure 3.1: Illustration of MPC with Contingency Predictions
Remark 3.9. The tightening of the constraint sets in (3.34)-(3.39) is essentially carried out twice,
to account for the effect of the disturbance before and after contingency activation.
31
3. FEASIBLE CONTINGENCIES
Algorithm 3.1: MPC with Feasible Contingency
1 while x(k) /∈ T do
2 Solve (3.40) subject to (3.41);
3 Apply the first element of the resulting optimal input sequence to the system;
4 k← k + 1;
5 if contingency α requires activation then
6 Switch to an auxiliary controller to handle the contingency;
7 exit;
8 end
9 end
Remark 3.10. The states and inputs after anticipated contingency activation are only con-
strained for j ≥ 1 as the current state, which corresponds to j = 0, cannot be modified within
the MPC optimisation problem. However, for k > 0, a feasible trajectory for j = 0 will exist, as
will be subsequently shown.
Algorithm 3.1 describes how this optimisation problem is used to control the system to the
target set, whilst ensuring that contingency α is available if required.
Theorem 3.11 (Robust Recursive Feasibilty). Given the optimal solution to the MPC optimisa-
tion in Algorithm 3.1 at time k, namely the sequences {u∗(k + j|k)}, {x∗(k + j|k)}, {y∗(k + j|k)},
{φ∗(α)j (i|k)}, {µ∗(α)j (i|k)} and {ψ∗(α)j (i|k)}, as well as the completion times N∗(k) and M∗(α)j (k), a
feasible solution to the optimisation problem at time k + 1 is given by
uˆ(k + j + 1|k + 1) = u∗(k + j + 1|k) + P(j)w(k), ∀j ∈ Z[0,N∗(k)−2] (3.43a)
xˆ(k + j + 1|k + 1) = x∗(k + j + 1|k) + L(j)w(k), ∀j ∈ Z[0,N∗(k)−1] (3.43b)
yˆ(k + j + 1|k + 1) = y∗(k + j + 1|k) + Q(j)w(k), ∀j ∈ Z[0,N∗(k)−2] (3.43c)
µˆ
(α)
j (i|k + 1) = µ∗(α)j+1 (i|k) + P˜(α)(i)L(j)w(k), ∀j ∈ Z[0,N∗(k)−2], i ∈ Z[0,M∗(α)j+1 (k)−1] (3.43d)
φˆ
(α)
j (i|k + 1) = φ∗(α)j+1 (i|k) + L˜(α)(i)L(j)w(k), ∀j ∈ Z[0,N∗(k)−2], i ∈ Z[0,M∗(α)j+1 (k)] (3.43e)
ψˆ
(α)
j (i|k + 1) = ψ∗(α)j+1 (i|k) + Q˜(α)(i)L(j)w(k), ∀j ∈ Z[0,N∗(k)−2], i ∈ Z[0,M∗(α)j+1 (k)−1] (3.43f)
Nˆ(k) = N∗(k)− 1 (3.43g)
Mˆ(α)j (k + 1) = M
∗(α)
j+1 (k). (3.43h)
Proof. The proof follows a similar methodology to [RH06a], albeit modified to include the
contingency prediction variables and disturbance-feedback constraint tightening.
Satisfaction of the initial conditions (3.22) and (3.25) can be proven as follows. For the primary
trajectory prediction, (3.23) and Algorithm 3.1 can be used to show that
xˆ(k + 1|k + 1) = x∗(k + 1|k) + w(k) (3.44)
32
3.3. ROBUST MPC WITH CONTINGENCIES
= Ax(k) + Bu∗(k|k) + w(k) (3.45)
= Ax(k) + Bu(k) + w(k) (3.46)
= x(k + 1). (3.47)
In a similar way, (3.43b) and (3.43e) can be used to show that
φˆ
(α)
j (0|k + 1) = φ∗(α)j+1 (0|k) + L˜(α)(0)L(j)w(k) (3.48)
= x∗(k + j + 1|k) + L(j)w(k) (3.49)
= xˆ(k + j + 1|k + 1). (3.50)
For the system dynamics, (3.23) and (3.29) can be used to show that the candidate solution
satisfies
xˆ(k + j + 2|k + 1) = x∗(k + j + 2|k) + L(j + 1)w(k) (3.51)
= Ax∗(k + j + 1|k) + Bu∗(k + j + 1|k) + L(j + 1)w(k) (3.52)
= A(xˆ(k + j + 1|k + 1)− L(j)w(k)) + B(uˆ(k + j + 1|k + 1)− P(j)w(k))
+ L(j + 1)w(k)
(3.53)
= Axˆ(k + j + 1|k + 1) + Buˆ(k + j + 1|k + 1)− (AL(j) + BP(j))w(k)
+ L(j + 1)w(k)
(3.54)
= Axˆ(k + j + 1|k + 1) + Buˆ(k + j + 1|k + 1)− L(j + 1)w(k)
+ L(j + 1)w(k)
(3.55)
= Axˆ(k + j + 1|k + 1) + Buˆ(k + j + 1|k + 1). (3.56)
For the contingency predictions, (3.26) and (3.32) show that
φˆ
(α)
j (i + 1|k + 1) = φ∗(α)j+1 (i + 1|k) + L˜(α)(i + 1)L(j)w(k) (3.57)
= A˜(α)φ∗(α)j+1 (i|k) + B˜(α)µ∗(α)j+1 (i|k) + L˜(α)(i + 1)L(j)w(k) (3.58)
= A˜(α)(φˆ(α)j (i|k + 1)− L˜(α)(i)L(j)w(k)) + B˜(α)(µˆ(α)j (i|k + 1)
− P˜(α)(i)L(j)w(k)) + L˜(α)(i + 1)L(j)w(k)
(3.59)
= A˜(α)φˆ(α)j (i|k + 1) + B˜(α)µˆ(α)j (i|k + 1)
− (A˜(α) L˜(α)(i) + B˜(α)P˜(α)(i))L(j)w(k) + L˜(α)(i + 1)L(j)w(k)
(3.60)
= A˜(α)φˆ(α)j (i|k + 1) + B˜(α)µˆ(α)j (i|k + 1)− L˜(α)(i + 1)L(j)w(k)
+ L˜(α)(i + 1)L(j)w(k)
(3.61)
= A˜(α)φˆ(α)j (i|k + 1) + B˜(α)µˆ(α)j (i|k + 1). (3.62)
For the output dynamics, (3.24) and (3.30) can be used to show that the candidate solution
33
3. FEASIBLE CONTINGENCIES
satisfies
yˆ(k + j + 1|k + 1) = y∗(k + j + 1|k) + Q(j)w(k) (3.63)
= Cx∗(k + j + 1|k) + Du∗(k + j + 1|k) + Q(j)w(k) (3.64)
= C(xˆ(k + j + 1|k + 1)− L(j)w(k)) + D(uˆ(k + j + 1|k + 1)− P(j)w(k))
+ Q(j)w(k)
(3.65)
= Cxˆ(k + j + 1|k + 1) + Duˆ(k + j + 1|k + 1)− (CL(j) + DP(j))w(k)
+ Q(j)w(k)
(3.66)
= Cxˆ(k + j + 1|k + 1) + Duˆ(k + j + 1|k + 1)−Q(j)w(k) + Q(j)w(k)
(3.67)
= Cxˆ(k + j + 1|k + 1) + Duˆ(k + j + 1|k + 1) (3.68)
For the contingency predictions, (3.27) and (3.33) show that
ψˆ
(α)
j (i|k + 1) = ψ∗(α)j+1 (i|k) + Q˜(α)(i)L(j)w(k) (3.69)
= C˜(α)φ∗(α)j+1 (i|k) + D˜(α)µ∗(α)j+1 (i|k) + Q˜(α)(i)L(j)w(k) (3.70)
= C˜(α)(φˆ(α)j (i|k + 1)− L˜(α)(i)L(j)w(k)) + D˜(α)(µˆ(α)j (i|k + 1)
− P˜(α)(i)L(j)w(k)) + Q˜(α)(i)L(j)w(k)
(3.71)
= C˜(α)φˆ(α)j (i|k + 1) + D˜(α)µˆ(α)j (i|k + 1)− (C˜(α) L˜(α)(i) + D˜(α)P˜(α)(i))L(j)w(k)
+ Q˜(α)(i)L(j)w(k)
(3.72)
= C˜(α)φˆ(α)j (i|k + 1) + D˜(α)µˆ(α)j (i|k + 1)− Q˜(α)(i)L(j)w(k) + Q˜(α)(i)L(j)w(k)
(3.73)
= C˜(α)φˆ(α)j (i|k + 1) + D˜(α)µˆ(α)j (i|k + 1) (3.74)
Satisfaction of the output constraints can be verified by using (3.38) and the definition of the
P-difference (2.55).
y∗(k + j + 1|k) ∈ Y(j + 1) = Y(j)	Q(j)W
=⇒ y∗(k + j + 1|k) + Q(j)W = yˆ(k + j + 1|k + 1) ∈ Y(j). (3.75)
Similarly,
ψ
∗(α)
j+1 (i|k) ∈ Y˜ (α)j+1(i) = Y˜ (α)j (i)	 Q˜(α)(i)L(j)W
=⇒ ψ∗(α)j+1 (i|k) + Q˜(α)(i)L(j)w(k) = ψˆ(α)j (i|k + 1) ∈ Y˜ (α)j (i). (3.76)
34
3.3. ROBUST MPC WITH CONTINGENCIES
In the same manner, (3.39) can be used to verify that the candidate solution satisfies the terminal
constraints
x∗(k + N∗(k)|k) ∈ T (N∗(k)) = T (N∗(k)− 1)	 L(N∗(k)− 1)W
=⇒ x∗(k + N∗(k)|k) + L(N∗(k)− 1)W = xˆ(k + Nˆ(k + 1)|k + 1) ∈ T (Nˆ(k + 1)). (3.77)
Similarly
φ
∗(α)
j+1 (M
∗(α)
j+1 (k)|k) ∈ T˜ (α)j+1 (M∗(α)j+1 (k)) = T˜ (α)j (M∗(α)j+1 (k))	 L˜(α)(M∗(α)j+1 (k))L(j)W
=⇒ φ∗(α)j+1 (M∗(α)j+1 (k)|k) + L˜(α)(M∗(α)j+1 (k))L(j)w(k)
= φˆ
(α)
j (Mˆ
(α)
j (k + 1)|k + 1) ∈ T˜ (α)j (Mˆ(α)j (k + 1)). (3.78)

By induction, this means that, given an initial feasible solution at time k, there exists a feasible
solution until completion is reached, in the absence of contingency activation. The case of
contingency activation is considered in the following theorem.
Theorem 3.12 (Guaranteed Contingency Feasibility). Given an initial feasible solution, Algorithm
3.1 ensures that, if contingency α is activated at time ks, there exists a control law that drives the state
of the system to T˜ (α) in M¯ steps or fewer, for any admissible sequence of disturbances after contingency
activation.
Proof. Given an initial feasible solution, Theorem 3.11 guarantees that there exists a nominal
input sequence
µ
∗(α)
0 (i|ks), i ∈ Z[0,M∗(α)0 (ks)−1] (3.79)
and resulting trajectory to the contingency set at any time 0 < ks < kc that satisfies
ψ
∗(α)
0 (i|ks) ∈ Y˜ (α)0 (i) (3.80)
φ
∗(α)
0 (M
∗(α)
0 (ks)|ks) ∈ T˜ (α)0 (M∗(α)0 (ks)), (3.81)
where M∗(α)0 (ks) ≤ M¯(α) as enforced by (3.41). Notice also that the sets Y˜ (α)0 (i) and T˜ (α)0 (i)
have been defined to ensure that, given the feasible input at time ks, a feasible input to the
system at time ks + i that guarantees satisfaction of the original (i.e. non-tightened) output
constraints ψ(i) ∈ Y˜ and terminal constraint φ(M∗(α)0 (ks)) ∈ T˜ (α) is given by
µ(i) =
µ
∗(α)
0 (0|ks) i = 0
µ
∗(α)
0 (i|ks) +∑i−1l=0 P˜(α)(i− l − 1)ω(l) 0 < i < M(α)0 (ks).
(3.82)
This result follows from modifying Theorem 3.11 to consider the contingency trajectory, and
applying it inductively to show to show that there exists a feasible solution at time ks + i + 1
given any feasible (not necessarily optimal) solution at time ks + i.
35
3. FEASIBLE CONTINGENCIES
Using q to index the predictions that would be made after contingency activation, the predicted
optimal contingency input trajectory at time ks + i is denoted {µ∗(α)0 (i + q|ks)}. Then, a feasible
predicted input trajectory at time ks + 1 is given by
uˆ(ks + 1+ q|ks + 1) = µ∗(α)0 (q + 1|ks) + P˜(α)(q)ω(0), q ∈ Z[0,M∗(α)0 (ks)−2] (3.83)
Now, assume for any i that a feasible input trajectory is given by
uˆ(ks + i + q|ks + i) = µ∗(α)0 (q + i|ks) +
i−1
∑
l=0
P˜(α)(q + i− l − 1)ω(l), q ∈ Z
[0,M∗(α)0 (ks)−i−1]
(3.84)
Then applying the modified Theorem 3.11 to this trajectory shows that
uˆ(ks + i + 1+ q|ks + i + 1) = uˆ(ks + i + 1+ q|ks + i) + P˜(α)(q)ω(i) (3.85)
= µ
∗(α)
0 (q + 1+ i|ks) +
i−1
∑
l=0
P˜(α)(q + i− l)ω(l) + P˜(α)(q)ω(i)
(3.86)
= µ
∗(α)
0 (q + 1+ i|ks) +
i
∑
l=0
P˜(α)(q + i− l)ω(l). (3.87)
where q ∈ Z
[0,M∗(α)0 (ks)−i−1]
. Given that the MPC algorithm will, at each time step, apply the
first input, that is, at q = 0, the actual inputs are found by substituting this value into (3.84) for
each i.
Hence, (3.82) defines a time-varying affine control law that applies direct feedback to the
disturbances ω( · ) entering after contingency activation. This implies that there exists a control
law that drives the system from the current state to the target set in M¯ steps or fewer, for any
admissible disturbance realisation, as required by Criterion 3.7. 
Remark 3.13. Theorem 3.12 provides a candidate control law to be used under contingency
activation. This means that, in the event of an optimiser failure, the law could actually be
applied to the system to guide it to the contingency terminal set.
Theorem 3.14 (Robust Finite-Time Completion). If the cost weighting matrix Θ in (3.19) is chosen
to satisfy η > 0, where
η = 1− ‖Θ‖max
w∈W
N¯−2
∑
j=0
‖Q(j)w‖ , (3.88)
the target set will be reached in at most H steps, where
H =
⌊
J∗(x(0))
η
⌋
. (3.89)
Proof. Denote the cost of the candidate solution (3.43) as Jˆ( · ). Then, applying the triangle
36
3.3. ROBUST MPC WITH CONTINGENCIES
inequality and maximising over all allowable disturbances and horizon lengths,
J∗(x(k))− Jˆ(x(k + 1)) =
N∗(k)−1
∑
j=0
(1+ ‖Θ(y∗(k + j|k)− yr)‖)
−
N∗(k)−2
∑
j=0
(1+ ‖Θ(y∗(k + j + 1|k) + Q(j)w(k)− yr)‖)
(3.90)
≥ 1−
N∗(k)−2
∑
j=0
‖ΘQ(j)w(k)‖ (3.91)
≥ 1− ‖Θ‖max
w∈W
N¯−2
∑
j=0
‖Q(j)w‖ (3.92)
= η (3.93)
From optimality, the cost Jˆ( · ) is an upper bound on the optimal cost at time k + 1, hence
J∗(x(k))− J∗(x(k + 1)) ≥ J∗(x(k))− Jˆ(x(k + 1)) (3.94)
≥ η. (3.95)
The value function is therefore a Lyapunov-like function, which reduces at each time step by at
least η. As the cost is non-negative by construction, even within the terminal set, this implies
that completion must be reached in at most bJ∗(x(0))/ηc steps. 
Remark 3.15. Note that η > 0 is sufficient to guarantee finite-time completion, but is not a
necessary condition, as there may be other choices of cost weighting that still maintain the
guarantee for a given system.
With an appropriate choice of cost function after contingency activation, or simply applying
the control law detailed in Theorem 3.12, finite-time entry of the contingency target set is also
guaranteed.
3.3.1 MIP implementation
As stated in Chapter 2, the cost function in (3.40) is, in general, a non-convex function of
the horizon length. Mixed-integer programming is one method that can be used to solve
the non-convex optimisation problem and incorporate the horizon lengths before and after
switching as decision variables. The prediction model is augmented with the binary variable
sequences {b(k + j|k) ∈ {0, 1}}, {β(α)j (i|k) ∈ {0, 1}} and {v(α, j) ∈ {0, 1}}. The variables b( · )
and β(α)j ( · ) take the value 1 at the time at which completion is achieved on the primary and
contingency trajectory predictions respectively. The variables v( · , · ) indicate whether or not
the contingency α is enabled, (i.e. required to be feasible) at prediction step j, and will be used
in subsequent sections of this chapter.
37
3. FEASIBLE CONTINGENCIES
The binary variables can be used to represent the value function (3.40), where the summation
is now over a fixed horizon length. This function takes the form
J∗(x(k)) = min
θb
N¯−1
∑
j=0
(jb(k + j|k) + ‖Θ(y(k + j|k)− yr‖) , (3.96)
where the optimisation is now with respect to the set of decision variables
θb =
⋃
i,j
{
u(k + j|k), µ(α)j (i|k), b(k + j|k), β(α)j (i|k)
}
. (3.97)
The optimisation is carried out subject to the mixed logic and dynamics constraints
b(k + j|k) = 1 =⇒ x(k + j|k) ∈ T (j) (3.98)
β
(α)
j (i|k) = 1 =⇒ φ(α)j (i|k) ∈ T˜ (α)j (i) (3.99)
j
∑
l=0
b(k + l|k) = 0 =⇒ y(k + j|k) ∈ Y(j) (3.100)
j
∑
l=0
b(k + l|k) +
i
∑
l=0
β
(α)
j (l|k) + (1− v(α, j)) = 0 =⇒ ψ(α)j (i|k) ∈ Y˜ (α)j (i) (3.101)
N¯
∑
l=0
b(k + l|k) = 1 (3.102)
M¯(α)
∑
l=0
β
(α)
j (l|k) +
j
∑
l=0
b(k + l|k)− v(α, j) ≥ 0, (3.103)
where j ∈ Z[0,N¯] for the primary trajectory variables, j ∈N[1,N¯] for the contingency variables
and i ∈ Z[0,M¯(α)]. Note that the condition (3.102) ensures that completion occurs somewhere
along the primary trajectory and (3.103) guarantees that if primary completion has not occurred
and the contingency is active at step j, completion must be reached somewhere along the
contingency prediction propagated from x(k + j|k). Also note that the output-based operating
costs are relaxed once contingency completion is achieved by virtue of (3.101). If the constraint
sets are polyhedral, a “big-M” formulation can be used to implement these logical implications
in a similar way to Section 2.7. By replacing the optimisation problem and constraints in
Algorithm 3.1 with (3.96) and (3.98) - (3.103) respectively, the resulting optimisation problem is
a MIP.
3.4 Contingency Scenarios
The algorithm outlined in the previous section guarantees that a single contingency α is feasible.
This section extends the basic case to other contingency scenarios, namely costed contingen-
cies, multiple contingencies, prioritised contingencies and state dependent contingencies. It
shows how these scenarios can be formulated using the mixed-integer framework previously
38
3.4. CONTINGENCY SCENARIOS
developed.
3.4.1 Costed Contingencies
It is possible that there could be some freedom in the choice of the nominal trajectory. For an
appropriately sized set F˜ (α)M¯ ( · ), the formulation presented in the previous section will ensure
only the existence of a feasible trajectory to the contingency set, whilst attempting to make
the cost of the primary trajectory as small as possible. However, if the cost of the contingency
trajectory is also important, this could affect the shape of the primary trajectory. The value
function can be modified in two different ways to account for this.
Worst-case contingency cost
It is possible to minimise a weighted sum of the cost to the primary target and the worst case
cost for contingency α, from each state prediction. This is achieved by the value function
J∗I (x(k)) = min
θ
(
N(k)−1
∑
j=0
(1+ ‖Θ(y(k + j|k)− yr)‖)
+ ρ max
j∈N[1,N¯−1]
M(α)j (k)−1
∑
i=0
(
1+
∥∥∥Θ˜(α) (ψ(α)j (i|k)− ψ(α)r )∥∥∥)
 , (3.104)
where Θ˜(α) is an appropriately sized cost weighting matrix, ρ ∈ R+ is the weighting on the
contingency cost and ψ(α)r is some contingency output reference. Assuming that 1-norms are
used, this can be written in mixed-integer form as
min
θb,s(k)
N¯−1
∑
j=0
(jb(k + j|k) + ‖Θ(y(k + j|k)− yr)‖1) + ρs(k), (3.105)
where the slack variable s(k) ≥ 0 satisfies the constraint
s(k) ≥
M¯(α)−1
∑
i=0
(
iβ(α)j (i|k) +
∥∥∥Θ˜(α) (ψ(α)j (i|k)− ψ(α)r )∥∥∥1) , ∀j ∈N[1,N¯−1]. (3.106)
The 1-norm in (3.106) can be represented by extra slack vectors to bound each component from
above and below, and subsequently summing these together. By optimality, this will provide a
tight bound.
Also by optimality, s(k) will be an upper bound only on those contingency trajectories that
originate prior to completion, since all output constraints are relaxed upon completion. It
can be shown that for an appropriate choice of cost weightings, the value function remains a
Lyapunov-like function, and thus preserves the guarantee of finite-time completion.
39
3. FEASIBLE CONTINGENCIES
Denoting the cost of the feasible solution (3.43) at time k + 1 by JˆI( · ), observe that
J∗I (x(k))− JˆI(x(k + 1))
= J∗(x(k))− Jˆ(x(k + 1)) + ρ max
j∈N[1,N∗(k)−1]
M∗(α)j (k)−1
∑
i=0
(
1+
∥∥∥Θ˜(α)(ψ∗(α)j (i|k)− ψ(α)r )∥∥∥)
− ρ max
j∈N[1,N∗(k)−2]
M∗(α)j+1 (k)−1
∑
i=0
(
1+
∥∥∥Θ˜(α) (ψ∗(α)j+1 (i|k) + Q˜(α)(i)L(j)w(k)− ψ(α)r )∥∥∥)
(3.107)
≥ η + ρ max
j∈N[1,N∗(k)−1]
M∗(α)j (k)−1
∑
i=0
(
1+
∥∥∥Θ˜(α) (ψ∗(α)j (i|k)− ψ(α)r )∥∥∥)
− ρ max
j∈N[1,N∗(k)−1]
M∗(α)j+1 (k)−1
∑
i=0
(
1+
∥∥∥Θ˜(α) (ψ∗(α)j+1 (i|k)− ψ(α)r )∥∥∥)
− ρ max
j∈N[1,N∗(k)−2]
M∗(α)j+1 (k)−1
∑
i=0
∥∥∥Θ˜(α)Q˜(α)(i)L(j)w(k)∥∥∥
(3.108)
≥ η − ρ
∥∥∥Θ˜(α)∥∥∥ max
j∈N[1,N¯−2]
w∈W
M¯(α)−1
∑
i=0
∥∥∥Q˜(α)(i)L(j)w∥∥∥ . (3.109)
Hence, by choosing the cost weightings Θ, Θ˜(α) and ρ such that the RHS of (3.109) is greater
than zero, the finite-time completion guarantee still holds.
Mean contingency cost
In addition to the primary trajectory cost, the mean contingency cost can also be penalised with
a value function of the form
J∗II(x(k)) = min
θ
(
N(k)−1
∑
j=0
(1+ ‖Θ(y(k + j|k)− yr)‖)
+
ρ
(N(k)− 1)
N(k)−1
∑
j=1
M(α)j (k)−1
∑
i=0
(
1+
∥∥∥Θ˜(α) (ψ(α)j (i|k)− ψ(α)r )∥∥∥)
 . (3.110)
With the addition of extra variables, this can be cast into mixed-integer form if a 1-norm cost
function is used. Define
σ1(k) ≥
M¯(α)−1
∑
i=0
(
iβ(α)1 (i|k) +
∥∥∥Θ˜(α)(ψ(α)1 (i|k)− ψ(α)r )∥∥∥1) (3.111)
σj(k) ≥ σj−1(k) + 1j
(
M¯(α)−1
∑
i=0
(
iβ(α)j (i|k) +
∥∥∥Θ˜(α) (ψ(α)j (i|k)− ψ(α)r )∥∥∥1)− σj−1(k)
)
, (3.112)
40
3.4. CONTINGENCY SCENARIOS
for j ∈N[2,N¯−1]. Then by using a value function of the form
min
θb,s(k)
N¯−1
∑
j=0
(jb(k + i|k) + ‖Θ(y(k + j|k)− yr)‖1) + ρs(k), (3.113)
where the slack variable s(k) ≥ 0 is redefined to satisfy the logical implications
b(k + j + 1|k) = 1 =⇒ s(k) ≥ σj(k), ∀j ∈N[1,N¯−1], (3.114)
the mean contingency cost is now penalised in addition to the normal operating cost. Slack
vectors can again be used to replace the 1-norm in (3.111)-(3.112). In this case the optimal cost
J∗II( · ) cannot be used as a Lyapunov-like function, as the mean is not, in general, non-increasing.
However, it is easy to show that
J∗II(x(k)) ≤ J∗I (x(k)), (3.115)
for any x(k). Hence, by choosing ρ and the cost weightings to satisfy the conditions (3.109)
for the worst-case contingency cost, finite time completion is also guaranteed for the mean
contingency cost, as an upper bound on J∗II( · ) reduces by at least some fixed amount at each
time step.
Remark 3.16. The cost weightings on the contingencies only consider nominal costs. Under
disturbances, the actual closed loop cost after contingency activation will differ. This is, of
course, true of the primary trajectory cost as well.
3.4.2 Multiple & Prioritised Contingencies
Algorithm 3.1 can be extended in a straightforward manner to ensure feasibility of multiple
contingencies. To do so, extra prediction variables need to be augmented to the system.
Specifically, sequences of states, inputs and outputs are required for each α ∈ N[1,nc]. Then
constraints will be enforced on each of these sequences according to (3.41). Define the set
C ⊆N[1,nc], indicating the indices of the contingencies that are required to be available. Then,
by enforcing the constraint
v(α, j) = 1, ∀α ∈ C, j ∈N[1,N¯−1], (3.116)
there must exist feasible trajectories from each primary state prediction in the event of any
of these contingencies being activated. This may, of course, make finding an initial feasible
solution more difficult.
Another manner in which the contingencies could be used is to ensure the availability of at
least one contingency within the set C at all times 0 < k < ks. By enforcing the constraints
∑
α∈C
v(α, j) ≥ 1, ∀j ∈N[1,N¯−1] (3.117)
41
3. FEASIBLE CONTINGENCIES
at least one of the contingencies in C will be available at each time step. Furthermore, by
modifying the cost function, it is also possible to prioritise the availability of contingencies
at a given prediction step if multiple options are feasible. Hence, a contingency that makes
the primary trajectory cost worse can be given a higher priority if demanded by operational
requirements. The new value function is given by
J∗III(x(k)) = min
θb,v( · , · )
N¯−1
∑
j=0
(
jb(k + i|k) + ‖Θ(y(k + j|k)− yr)‖+ ∑
α∈C
γαv(α, j)
)
, (3.118)
where γα ∈ R+ is the weighting on contingency α. Note that this preserves the guarantee
of finite-time completion, as the feasible solution at the next time step still has the same
contingency availability, so the cost term involving the variables v( · , · ) reduces.
3.4.3 State-Dependent Contingencies
It may be the case that a particular contingency is required to be available depending on the
current state of the system. Given the output constraint set Y , the output-admissible set of
states is given by
O = {x | ∃u : Cx + Du ∈ Y}. (3.119)
Now define a set A ⊆ O, which is a set that defines a region of the state space in which one
or more contingencies are required to be available. A partition of A is then specified by the
sequence of sets {Aα}α∈C , which satisfy⋃
α∈C
Aα = A (3.120)
Aα1 ∩Aα2 = ∅, ∀α1, α2 ∈ C, α1 6= α2. (3.121)
The ideal requirement is that contingency α is available if the state of the system is within Aα.
However, in order to maintain recursive feasibility, the sets Aα need to be enlarged. This is due
to the fact that if x∗(k + j+ 1|k) ∈ Aα0 , for some α0 there is no guarantee that xˆ(k + j+ 1|k + 1)
will also lie in Aα0 , due to the action of the disturbance at time k. Therefore, define the enlarged
sets
A¯α = Aα ⊕
N¯−1⊕
j=0
L(j)W . (3.122)
Then, recursively define the sets
A¯α(0) = A¯α (3.123)
A¯α(j + 1) = A¯α(j)	 L(j)W . (3.124)
By now enforcing the constraints
x(k + j|k) ∈ A¯α(j) =⇒ v(α, j) = 1, (3.125)
42
3.5. SIMPLIFICATIONS
the availability of the appropriate contingency in the MPC optimisation problem is ensured,
depending on the predicted state. This is due to the fact that, from the design of sets Aα
as a partition of the entire feasible state-space, and given that the relationship between the
Pontryagin difference and Minkowski sum [KG95] ensures that A¯α(j) ⊇ Aα for all j ∈ Z[0,N¯],
the system state is always guaranteed to be in at least one of the sets A¯α(j). This will, however,
require multiple contingencies to be available in regions of overlap.
Remark 3.17. The implication (3.125) is reversed from all of the others encountered so far. Extra
binary variables are needed to represent this reverse implication. If the sets can be represented
by r half space constraints as
A¯α(j) =
{
x
∣∣∣ G(α)(j)x ≤ h(α)(j), G(j) ∈ Rr×n, h(j) ∈ Rn×1} , (3.126)
by introducing the binary vector δ(α)(j) ∈ {0, 1}r×1, the implication can be expressed in
mixed-integer form as
G(α)(j)x ≥ f (α)(j)−M
(
1− δ(α)(j)
)
(3.127)
1− v(α, j) ≤
r
∑
i=1
δ
(α)
i (j), (3.128)
for all j ∈ ZN¯−10 , where δ(α)i (j) denotes the ith binary element of δ(α)(j) and M is a suitably
large number.
Recursive feasibility of the resulting optimisation problem is then readily apparent, as
x∗(k + j + 1|k) ∈ A¯α(j + 1) =⇒ x∗(k + j + 1|k) + L(j)W = xˆ(k + j + 1|k + 1) ∈ A¯α(j).
(3.129)
Hence, a shifted version of the optimal solution at time k with the disturbance feedback
term added still requires the same contingency availability. This essentially means that the
disturbance feedback term doesn’t cause the state to “jump” into another region where a
different contingency is required to be feasible, or worse still, move into the set A from being
outside it.
Remark 3.18. It is possible to combine costed contingencies with multiple contingency availab-
ility and state dependent contingencies, through appropriate addition of prediction variables
and modification of the cost function.
3.5 Simplifications
The formulation presented in the previous two sections uses additional prediction variables for
contingency dynamics, as well as integer variables for indicating completion on the primary
and contingency trajectories. This can become very computationally expensive if there are
many possible contingencies. There are three simplifications that can be applied to reduce
complexity.
43
3. FEASIBLE CONTINGENCIES
3.5.1 Fixed Contingency Horizon
By assuming a fixed horizon on the contingency trajectories, the need for binary variables on
these trajectory predictions is eliminated. However, this may make it more difficult to find an
initially feasible solution, and will produce a more conservative primary trajectory. A fixed
contingency horizon may be appropriate when the contingency terminal set T˜ (α)0 is robust
control invariant and upon contingency activation, the control problem is to remain within the
contingency terminal set or stabilise to a setpoint within this set.
3.5.2 Explicit Contingency Feasibility
If a fixed horizon is used, it is possible to further simplify the problem by explicitly calculating
the set K˜(α)M (T˜ (α)), assuming that the target set after contingency activation is to be entered in
M(α) steps. For the variable horizon case, it is possible to calculate the set F˜ (α)M (T˜ (α)), but this
will, in general, be non convex, as it will be the union of the controllable sets for each possible
horizon length after contingency activation. As each state on the primary trajectory only has to
be in one of these sets, binary variables will be needed to enforce the constraints.
It may be asked why the online solution is ever required if the robust controllable sets for
contingency feasibility can be computed. After all, even with explicit variable horizon con-
tingencies, the resulting mixed-integer program contains fewer decision variables. However,
without parameterising the contingency trajectories, it is not possible to place costs on them
and incorporate these into the optimisation problem. Furthermore, with non-convex collision
avoidance constraints, many sets will need to be calculated for every combination of avoidance
constraints over the whole prediction horizon. Finally, the contingency predictions provide a
feasible solution from Theorem 3.12 that can be applied in the case of optimiser failure under
contingency activation.
3.5.3 Prediction Variable Reduction
When considering the availability of at least one contingency as enforced by the condition
(3.117), it may be the case that all contingencies in C have the same dynamics and operating
constraints, but different terminal sets. In this case, only a single set of prediction variables is
needed for the states, inputs and outputs under contingency activation. The binary variables
indicating completion cannot be reduced, however. They are still required to distinguish
completion amongst all of the contingency target sets.
3.6 Examples
Examples illustrating the different contingency scenarios are now presented. A simple double-
integrator system is introduced first, with a contingency requirement to provide fault-tolerance.
This system is used to demonstrate the effect of costed, prioritised and state-based contingencies.
44
3.6. EXAMPLES
A cart-pendulum is introduced next, with an emergency stop contingency. It demonstrates the
feasible contingency concept applied to a setpoint regulation problem, with a fixed contingency
horizon.
3.6.1 Double-Integrator Vehicle
Consider a unit point mass moving in two dimensions, actuated by orthogonal forces in
the horizontal and vertical directions. Treating the mass as a fictitious vehicle, fuel usage
is quantified as the integral of the absolute actuator forces over time. Hence, there are two
inputs representing the force components and five states representing the position, velocity
and fuel. In order to measure the fuel usage, the inputs are split to capture the positive and
negative components separately, giving four inputs in the resulting system, constrained to be
non-negative.
Discretising with a sampling frequency of 1 Hz gives the state-space matrices
A =

1 0 1 0 0
0 1 0 1 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
 , B =

0.5 −0.5 0 0
0 0 0.5 −0.5
1 −1 0 0
0 0 1 −1
−1 −1 −1 −1
 (3.130)
The C and D matrices are chosen to output all states and feed through all inputs of the system,
that is, such that y = (x, u). In addition to the positivity constraints on each individual input,
in order to restrict the maximum actuation, the inputs satisfy the operating constraints
u(k) ∈ U = {u | |u| < 1}, (3.131)
where | · | indicates the row-wise absolute value. The states satisfy the constraints
x(k) ∈ X =
{
x
∣∣∣∣ [−20 −20 −2.5 −2.5 0]T ≤ x ≤ [20 20 2.5 2.5 10]T} , (3.132)
which restrict the maximum speed in each direction to 2.5 and ensure that the fuel level is
always non-negative. The remaining constraints on position and an upper bound on fuel use
simply ensure that the set X is a polytope. The system is assumed to start with 8 units of fuel,
which limits the overall use of control authority (or impulse applied to the vehicle).
Using VH-MPC, define the value function as
J∗(x(k)) = min
u(k),N(k)
N(k)−1
∑
j=0
1+ 0.1 ‖u(k + j|k)‖1 . (3.133)
This function can be expressed in the form (3.40) by an appropriate choice of Θ and yr. The
45
3. FEASIBLE CONTINGENCIES
primary target is now defined by
T =
{
x
∣∣∣∣ [−1 8 −2.5 −2.5 0]T ≤ x ≤ [1 10 2.5 2.5 10]T} . (3.134)
This set restricts the position of the vehicle to a translated square of side length 2 and ensures
that the velocity components in each coordinate direction have a maximum magnitude of 2.5.
It also ensures that the fuel is non-negative upon entry. The upper bound on the fuel is again
to ensure that T is a polytope.
The system is known to be vulnerable to a certain fault, where it begins to leak fuel at the rate
of 10% of its current level at each time step. At the same time, it requires 10% more fuel to
produce an equivalent force in the y actuator. This gives the new system dynamics defined by
matrices
A˜ =

1 0 1 0 0
0 1 0 1 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 0.9
 , B˜ =

0.5 −0.5 0 0
0 0 0.5 −0.5
1 −1 0 0
0 0 1 −1
−1 −1 −1.1 −1.1
 . (3.135)
When this fault occurs, the system is required, within 4 time steps, to reach another target set,
subject to the same input and state constraints as under normal operation. The contingency set
is given by
T˜ (1) =
{
φ
∣∣∣∣ [4 3.5 −2.5 −2.5 0]T ≤ φ ≤ [6 5.5 2.5 2.5 10]T} , (3.136)
which is another square, subject to the same velocity and fuel constraints as during normal
operation. To test the efficacy of robust contingency availability, or robust fault-tolerance, a
wind disturbance force, defined by the set
W =
{
w
∣∣∣∣ |w| ≤ [0 0 0.2 0.2 0]T} , (3.137)
acts on the system both during normal operation and under fault.
Robustness is ensured by tightening these constraint sets according to (3.34) - (3.37). The process
is simplified by choosing nilpotent disturbance feedback policies for constraint tightening,
which greatly reduces the number of sets that need to be calculated. It can be verified that, for
the system operating normally (3.130), the policy
P(j) =

[ −0.5 0 −0.75 0 0.25
0.5 0 0.75 0 0.25
0 −0.5 0 −0.75 0.25
0 0.5 0 0.75 0.25
]
j = 0[ −0.5 0 −0.75 0 0
0.5 0 0.75 0 0
0 −0.5 0 −0.75 0
0 0.5 0 0.75 0
]
j = 1
04×5 j > 1,
renders the system degree-2 nilpotent, meaning that A2 + ABP(0) + BP(1) = 0. For the system
46
3.6. EXAMPLES
under fault, the policy
P˜(1)(i) =

[ −0.5 0 −0.75 0 0.2036
0.5 0 0.75 0 0.2036
0 −0.5 0 −0.75 0.2240
0 0.5 0 0.75 0.2240
]
i = 0[ −0.5 0 −0.75 0 0
0.5 0 0.75 0 0
0 −0.5 0 −0.75 0
0 0.5 0 0.75 0
]
i = 1
04×5 i > 1,
also gives degree-2 nilpotency.
Simulations of the controlled, discretised system are performed over the same 20 random
sequences of allowable disturbances, for the various contingency scenarios detailed below.
Single Contingency
Figure 3.2 shows the trajectories under normal operation and anticipated contingency activation
at each time step. The cost function used once the contingency is activated takes the same
form as (3.133), but has a sufficiently small cost weighting on the inputs (0.05) to ensure that
completion does occur within four steps. On the contingency trajectories, the ‘+’ symbols
indicate the time steps.
It is clear that the distribution of trajectories entering the primary target is skewed towards the
right, illustrating the effect of requiring robust reachability of T˜ (1) from any state prediction.
Figures 3.3 and 3.4 show the speeds and actuator forces respectively. These quantities satisfy
the imposed constraints, indicated by the dashed lines, in the presence of the disturbance force.
The fuel use is displayed in Figure 3.5.
Costed Contingency
Using the same system and contingency requirement, a weighting can be placed on the mean
contingency cost as in (3.110). Defining the contingency cost as above, the mean cost over all
contingency predictions is weighted by 1.5 and added to the primary trajectory cost. This
essentially gives preference to contingency trajectories that complete sooner than the four-step
maximum. The effect of applying this weighting is evident in Figure 3.6. It can be seen that the
primary trajectories have now moved closer to the contingency set.
Multiple Prioritised Contingencies
Define a second contingency target set as
T˜ (2) =
{
φ
∣∣∣∣ [4 6.5 −2.5 −2.5 0]T ≤ φ ≤ [6 8.5 2.5 2.5 10]T} . (3.138)
The system is now required to reach either T˜ (1) or T˜ (2) if the fault occurs. Applying no priorities
on the two contingency target sets gives the trajectories shown in Figure 3.7. When the original
47
3. FEASIBLE CONTINGENCIES
−2 0 2 4 6 8
0
1
2
3
4
5
6
7
8
9
10
T
T˜ (1)
x
y
 
 
Primary
Contingency
Figure 3.2: Trajectories with contingency
0 1 2 3 4 5 6 7
−2
−1
0
1
2
k
Ve
lo
ci
ty
 x
0 1 2 3 4 5 6 7
−2
−1
0
1
2
k
Ve
lo
ci
ty
 y
Figure 3.3: Primary trajectory velocities
48
3.6. EXAMPLES
0 1 2 3 4 5 6 7
−1
−0.5
0
0.5
1
k
In
pu
t F
or
ce
 x
0 1 2 3 4 5 6 7
−1
−0.5
0
0.5
1
k
In
pu
t F
or
ce
 y
Figure 3.4: Primary trajectory control inputs
0 1 2 3 4 5 6 7
4.5
5
5.5
6
6.5
7
7.5
8
k
Fu
el
Figure 3.5: Primary trajectory fuel use
49
3. FEASIBLE CONTINGENCIES
−2 0 2 4 6 8
0
1
2
3
4
5
6
7
8
9
10
T
T˜ (1)
x
y
 
 
Primary
Contingency
Figure 3.6: Trajectories with mean-costed contingency
contingency has a smaller cost weighting placed on it compared to the new contingency using
(3.118), contingency completion in this set is now preferred. This is illustrated in Figure 3.8.
State-dependent Contingencies
To clearly demonstrate the effect of a state-dependent contingency, define a new contingency
target set as
T˜ (3) =
{
φ
∣∣∣∣ [8 3.5 −2.5 −2.5 0]T ≤ φ ≤ [10 5.5 2.5 2.5 10]T} . (3.139)
Assuming the same fault dynamics, the system is now required to reach T˜ (3) in the event of
the fault occurring, but only if the system state lies within the set
A =
{
x ∈ O |
[
0 1 0 0
]
x ≤ 4
}
, (3.140)
which is the set of all output-admissible states such that the the y coordinate is less than 4.
After expanding this set appropriately and applying the extra constraints (3.125), the system is
once again simulated over 20 admissible disturbance realisations. As shown in Figure 3.9, the
trajectories initially deviate to the right, to guarantee that T˜ (3) can be robustly reached. Once
this requirement is lifted, they moves back towards the left. Figure 3.10 shows how fuel use is
restricted whilst the contingency is required to be available.
50
3.6. EXAMPLES
−2 0 2 4 6 8
0
1
2
3
4
5
6
7
8
9
10
T
T˜ (1)
T˜ (2)
x
y
 
 
Primary
Contingency
Figure 3.7: Trajectories for two contingency sets
−2 0 2 4 6 8
0
1
2
3
4
5
6
7
8
9
10
T
T˜ (1)
T˜ (2)
x
y
 
 
Primary
Contingency
Figure 3.8: Trajectories for two contingency sets with priority
51
3. FEASIBLE CONTINGENCIES
0 2 4 6 8 10
0
1
2
3
4
5
6
7
8
9
10
T
T˜ (3)
x
y
 
 
Primary
Contingency
Figure 3.9: Trajectories for state-dependent contingency
0 1 2 3 4 5 6 7
2
3
4
5
6
7
8
k
Fu
el
Figure 3.10: Fuel use for state-dependent contingency
52
3.6. EXAMPLES
3.6.2 Cart-Pendulum
To demonstrate the effect of contingency behaviour on a setpoint regulation problem, consider
the idealised cart-pendulum shown in Figure 3.11. Given the applied force u, cart mass M,
m
l
u
θ
r
M
Figure 3.11: Diagram of cart-pendulum system
pendulum mass m, pendulum length l, the cart position r and pendulum angle θ, the system
satisfies the nonlinear differential equations,
(M + m sin2 θ)r¨ + mlθ˙2 sin θ −mg sin θ cos θ = u (3.141)
((M + m)l −ml cos2 θ)θ¨ + mlθ˙2 sin θ cos θ − (M + m)g sin θ = u cos θ. (3.142)
Linearising about θ = 0 and converting to state-space form gives the system
r˙
θ˙
r¨
θ¨
 =

0 0 1 0
0 0 0 1
0 mM g 0 0
0 M+mMl g 0 0


r
θ
r˙
θ˙
+

0
0
1
M
1
Ml
 u. (3.143)
Taking M = m = 1 kg, l = 1 m and g = 9.81 m s−2, this linearised system is discretised with a
sampling period Ts = 0.2 s, producing the discrete-time system with state space matrices,
A =

1 0.2094 0.2 0.0136
0 1.4187 0 0.2272
0 2.2289 1 0.2094
0 4.4578 0 1.4187
 , B =

0.0207
0.0213
0.2136
0.2272
 , (3.144)
53
3. FEASIBLE CONTINGENCIES
having linearised state vector x(k) =
[
r(kTs) θ(kTs) r˙(kTs) θ˙(kTs)
]T
and input vector
u(k) = u(kTs)).
An LQR controller is designed to stabilise the pendulum about (r, θ) = (rr, 0), where rr is some
reference position. The gain matrix is given by
K f =
[
2.0245 −41.3743 3.1303 −10.6857
]
. (3.145)
The system is subject to the pointwise-in-time state constraints (with respect to the reference)
x(k)− xr ∈ X =
{
x
∣∣∣∣ [−10 −0.8 −0.5 −0.5]T ≤ x ≤ [10 0.8 0.5 0.5]T} , (3.146)
where xr =
[
rr 0 0 0
]T
, and input constraints
u(k) ∈ U = {u | ‖u‖∞ ≤ 4.0}. (3.147)
These constraints are tightened by a small margin defined by the set
W =
{
w
∣∣∣∣ |w| ≤ [0.001 0.001 0.0001 0.0001]T} (3.148)
to help mitigate the effect of linearisation error. Taking a 15-step prediction horizon, a 15-step
nilpotent policy is used for tightening constraints. Given that this policy is nilpotent, the
control objective is to steer the system state to the setR, which is the largest nominally control
invariant set with respect to K f inside the terminal set given by
T =
{
x
∣∣∣∣ |x− xr| ≤ [0.1 0.5 0.5 0.5]T} . (3.149)
That is, for all x ∈ R ⊆ T ,
(A + BK f )(x− xr) + xr ∈ R (3.150)
K f (x− xr) ∈ U . (3.151)
Then, a VH-MPC controller will be used for transient control and once the terminal set is
entered, the controller K f will be applied, which is a dual mode MPC scheme. The value
function takes the form
J∗(x(k)) = min
u(k),N(k)
N(k)−1
∑
j=0
1+ 0.01 ‖u(k)‖1 . (3.152)
During the transition to the terminal set, a three-step fixed horizon emergency stop contingency
is required at each sampling instant. This contingency must keep the pendulum balanced
at the end of the emergency-stop manœuvre. Using the same disturbance set and constraint
tightening policy, the contingency requires the pendulum state to be steered into the largest
54
3.7. CONCLUSION
robust control invariant set R˜(1) with respect to K f and the disturbance L(2)W that lies within
the contingency terminal set
T (1) =
{
φ
∣∣∣∣ |φ− φs| ≤ [10 0.1 0.05 0.05]T} , (3.153)
where φs =
[
rs 0 0 0
]T
is any stopping reference state and the stopping position rs satisfies
|rs − rr| ≤ 10. This allows the cart at the time of contingency completion to be arbitrarily
positioned within a specified range of stopping positions relative to the desired cart position.
The contingency invariant set therefore ensures that, for all φ ∈ R˜(1) ⊆ T (1), there exists some
rs such that
(A + BK f )(φ− φs) + L(2)W + φs ∈ R˜(1) (3.154)
K f (φ− φs) ∈ U . (3.155)
This set places a more stringent terminal requirement on the pendulum angle and angular
velocity, as might be expected with an emergency stop constraint. Then, starting with the zero
state, a step input is used as the reference (i.e. rr = 1). The state trajectories for the nonlinear
pendulum dynamics with and without the imposition of contingency availability are shown in
Figures 3.12 and 3.13, for the cart and pendulum respectively. The control inputs are displayed
in Figure 3.14. Zero order hold discretisation is applied to the continuous states which are then
fed into the MPC controller.
Remark 3.19. The set R˜(1) ensures that there exists some rs at which the cart can be stabilised.
To find a suitable value for this position upon contingency activation, a quadratic program can
be solved to determine the φs closest to the terminal predicted contingency state at the current
time.
It can be seen that the contingency requirement produces a vastly different closed-loop be-
haviour. The only signal that approaches a constraint boundary is the cart velocity without
contingency, and this constraint is respected in spite of the nonlinearity as a result of to the
tightening applied. It can be seen that when the contingency is required to be available, the
cart velocity, pendulum angular velocity and input force are all reduced.
3.7 Conclusion
This chapter has introduced the notion of feasible contingencies and specified the criterion for a
contingency to be robustly feasible. It has demonstrated how a variable horizon controller can
be designed for linear systems to guarantee pointwise-in-time robust contingency availability.
This controller has been shown to preserve the guarantees of robust recursive feasibility and
robust finite-time completion with an appropriate choice of cost weightings. The basic formu-
lation has then been extended to show how these guarantees can be maintained with costed,
55
3. FEASIBLE CONTINGENCIES
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0
0.2
0.4
0.6
0.8
1
t
r
 
 
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
−0.2
0
0.2
0.4
0.6
0.8
t
dr
/d
t
Without Contingency
With Contingency
Setpoint
Figure 3.12: Step response for cart states
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
−0.1
−0.05
0
0.05
t
θ
 
 
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
−0.3
−0.2
−0.1
0
0.1
0.2
t
dθ
/d
t
Without Contingency
With Contingency
Figure 3.13: Step response for pendulum states
56
3.8. ACKNOWLEDGEMENTS
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
−4
−3
−2
−1
0
1
2
3
4
t
u
 
 
Without Contingency
With Contingency
Figure 3.14: Control input to cart
multiple, prioritised and state-dependent contingencies. Two example systems demonstrate
the power of the formulation.
Future work will consider the issue of further contingency activations along the contingency
trajectories, as well as techniques for reducing the computational complexity of the approach,
which depends on the use of mixed-integer programming in implementation. For the latter,
this could include the use of move blocking strategies described in Chapter 4, or intelligently
reducing the predicted states at which contingency availability is enforced.
3.8 Acknowledgements
This chapter builds upon a paper presented by the author at the 18th IFAC World Congress
[SM11]. The simulations in this chapter were carried out using MATLAB R2011b, YALMIP
[Löf04], CPLEX [CPL09] and the Invariant Set Toolbox [Ker00].
57
58
Chapter 4
Move Blocking for Complexity
Reduction
4.1 Introduction
A major issue with Variable Horizon MPC, as mentioned in Chapter 2, is that the cost function
optimised at each time step is not, in general, a convex function of the horizon length. Recall
that there are two common strategies for achieving global optimality in spite of this non-
convexity: evaluating multiple fixed horizon problems at each iteration or solving a mixed-
integer programming problem. With both methods, a bound must be placed on the maximum
horizon length to ensure a finite number of decision variables. However, this length must be
long enough to encompass the entire control manœuvre to completion, potentially requiring a
large number of decision variables for a long range control problem.
Move blocking [Mac02; Cag+07] presents a way of reducing the number of decision variables,
by curtailing the number of degrees of freedom in the MPC optimisation problem. It achieves
this by constraining groups of adjacent-in-time input decision variables to have the same value.
The manner in which such constraints are applied over time is termed the blocking regime.
In order to retain recursive feasibility guarantees, the blocking regime must be time-varying
[Cag+07], or additional constraints need to be enforced on the first prediction step [GI07;
GIK09].
In this chapter, a new formulation of VH-MPC is presented, utilising time-varying move
blocking to reduce complexity. It allows flexibility in the choice of blocking regime and,
through the use of tightened constraints, guarantees robust recursive feasibility and finite-time
completion in the presence of bounded disturbances. By using constraint tightening, the
formulation is, in general, computationally simpler than solving for blocked affine feedback
policies as outlined in [Old+09].
The contents of the chapter are structured as follows. Section 4.2 poses the control problem for
a linear discrete-time system. Section 4.3 introduces move blocking notation and various utility
59
4. MOVE BLOCKING FOR COMPLEXITY REDUCTION
functions. The controller formulation is described in Section 4.4. An example move blocking
regime and simulation results are presented in Section 4.6.
4.1.1 Nomenclature
Define the horizontal concatenation operator for matrices A(i) having the same row dimension
as
b
C
i=a
A(i) =
[
A(a) A(a + 1) · · · A(b)
]
, (4.1)
where a, b ∈ Z, a ≤ b. Also, define the direct sum of matrices C ∈ Rn×m and D ∈ Rp×q by
C⊕ D =
[
C 0n×q
0p×m D
]
. (4.2)
Note that this symbol has changed meaning from the previous chapter, where it represented
the Minkowski sum.
4.2 Problem Definition
Consider once again a discrete-time linear system with states x(k), inputs u(k) and outputs
y(k) specified by the difference equation
x(k + 1) = Ax(k) + Bu(k) (4.3)
y(k) = Cx(k) + Du(k), (4.4)
where A ∈ Rn×n, B ∈ Rn×m, C ∈ Rp×n and D ∈ Rp×m. It is assumed that the pair (A, B) is
controllable. The system is subject to pointwise-in-time input and state constraints of the form
u(k) ∈ U ⊂ Rm (4.5)
x(k) ∈ X ⊂ Rn, (4.6)
where the sets U and X are compact and closed respectively. The control objective is to steer the
state to a target set T ⊂ Rn, which is also closed, in time N < ∞, expressed as the constraint
x(N) ∈ T . (4.7)
The control input and horizon length are ideally chosen to minimise a p-norm cost function
that penalises inputs, states and time to completion. The value function takes the form
J∗(x(0)) = min
u,N
N−1
∑
k=0
(1+ ‖Θ(y(k)− yr)‖) , (4.8)
60
4.3. MOVE BLOCKING FORMULATION
for some matrix Θ ∈ Rp×p, sequence of inputs u and some reference output yr = Cxr + Dur,
where xr ∈ T and ur ∈ U respectively. The behaviour of the system after reaching T is not
specified as part of the control problem, and therefore there is no terminal cost in (4.8).
4.3 Move Blocking Formulation
In order to represent different move blocking regimes, blocking matrices will be employed.
These will be subsequently used in the development of a predictive controller to solve the
control problem posed in the previous section.
Definition 4.1. A blocking matrix S ∈ {0, 1}N×r is one that takes the form
S =

b1
b2
. . .
br
 =
r⊕
q=1
bq, (4.9)
where each bq = 1lq is termed a blocking vector, having length lq ∈N. The lengths satisfy
r
∑
q=1
lq = N. (4.10)
The set of all N × r blocking matrices will be denoted A(N, r).
Given a sequence of elements vi ∈ Rm×1, i ∈ N[1,N] and a corresponding sequence of block
values ωq ∈ Rm×1, q ∈ N[1,r], the blocking matrix S ∈ A(N, r) relates the two sequences
arranged in matrix form by the equation
N
C
i=1
vi =
(
r
C
q=1
ωq
)
ST =
(
C
q∈S
ωq
)
ST, (4.11)
where, with a slight abuse of notation, the symbol S has also been used to represent the ordered
set of all block indices, namelyN[1,r]. Now define the set-valued map of blocking matrices that
subdivide N elements or fewer into r blocks or fewer by
S(N, r) = {S ∈ A(i, q) | 0 ≤ i ≤ N, 0 ≤ q ≤ min{i, r}}. (4.12)
Definition 4.2. A relaxation S′ of a blocking matrix S ∈ A(N, r) is another blocking matrix in
which blocks can be amalgamated to give S. The set of all relaxations of S is given by
R(S) =
{
S′ ∈ RN×r′
∣∣∣ ∃M ∈ A(r′, r) : S = S′M, r′ ≥ r} . (4.13)
The identity blocking matrix (i.e. no blocking) is a relaxation of any blocking matrix. Each
blocking matrix is also a relaxation of itself. The set-valued mapR(S, r′) will be used to denote
61
4. MOVE BLOCKING FOR COMPLEXITY REDUCTION
all relaxations having exactly r′ blocks.
Define the function ` : A(N, r)×N[1,r] 7→N[1,N−r+1] for a given blocking matrix S ∈ A(N, r),
which returns the length of the qth block.
`(S, q) =
∥∥Seˆq∥∥1 , (4.14)
where eˆq is the qth Cartesian unit vector in r dimensions, (i.e. the qth column of Ir). Furthermore,
define a cumulative block length function Λ(S, q) as
Λ(S, q) =
0 q = 0∑qh=1 `(S, h), q ∈N[1,r]. (4.15)
A block-end indexing function is now defined in terms of Λ( · , · ).
µ(S, j) =

Λ(S, 1)− 1 Λ(S, 0) ≤ j < Λ(S, 1)
Λ(S, 2)− 1 Λ(S, 1) ≤ j < Λ(S, 2)
...
Λ(S, r)− 1 Λ(S, r− 1) ≤ j < Λ(S, r)
=
{
Λ(S, q)− 1 Λ(S, q− 1) ≤ j < Λ(S, q), q ∈N[1,r] . (4.16)
A function is also defined that returns the block number corresponding to the zero-based
element index.
σ(S, j) =

1 Λ(S, 0) ≤ j < Λ(S, 1)
2 Λ(S, 1) ≤ j < Λ(S, 2)
...
r Λ(S, r− 1) ≤ j < Λ(S, r)
=
{
q Λ(S, q− 1) ≤ j < Λ(S, q), q ∈N[1,r] . (4.17)
An example illustrating the action of functions Λ( · , · ), µ( · , · ) and σ( · , · ) can be found in
Figure 4.1.
Finally, define the shifted blocking matrix of S ∈ A(N, r), where N > 1 by the function
F(S) =

[
0(N−1)×1 IN−1
]
S `(S, 1) > 1[
0(N−1)×1 IN−1
]
S
01×(r−1)
Ir−1
 `(S, 1) = 1. (4.18)
Remark 4.3. Shifting essentially removes the first row of S and also the first column if `(S, 1) =
1 to ensure admissibility of the resulting matrix. This is a modification of the operation used in
[Cag+07] to suit the variable horizon.
62
4.3. MOVE BLOCKING FORMULATION
0 1 2 3 4 5 6 7 8
v1 v2 v3 v4 v5 v6 v7 v8 v9
1 2 3
2 6 9
1 1 5 5 5 5 8 8 8
j
vj+1
σ(S, j)
Λ(S, σ(S, j))
µ(S, j)
Figure 4.1: Values of functions Λ( · , · ) and µ( · , · ) for blocking matrix corresponding to N = 9
and block lengths {2, 4, 3}
Theorems 4.4 - 4.6 quantify the maximum number of blocking matrices and relaxations for a
given number of elements and blocks.
Theorem 4.4. Given N elements, the number of distinct blocking matrices comprising exactly r ≤ N
blocks is given by
|A(N, r)| =
(
N − 1
r− 1
)
. (4.19)
Proof. Any blocking matrix S ∈ A(N, r) maps bijectively onto a sequence of block lengths
l1, . . . , lr. These lengths are required to satisfy
r
∑
q=1
lq = N. (4.20)
Thus, a particular blocking matrix maps bijectively onto a particular composition [HM09] of N.
Using the notation
{l1, l2, . . . , lr} = {111 . . .1} (N terms), (4.21)
it can be seen that a unique composition can be obtained by replacing each ‘’ by either an
addition sign (+) or comma (,). This presents a binary choice, and since exactly r− 1 commas
must be placed to obtain r blocks, there are
( N−1
r−1
)
possible ways of doing this. 
Theorem 4.5. The number of blocking matrices that subdivide N or fewer elements into r or fewer
blocks is given by
|S(N, r)| = 2r−1 − 1+
N
∑
i=r
r
∑
q=1
(
i− 1
q− 1
)
. (4.22)
63
4. MOVE BLOCKING FOR COMPLEXITY REDUCTION
Proof. The result uses Theorem 4.4 and follows from expanding the summation
N
∑
i=1
min{i,r}
∑
q=1
(
i− 1
q− 1
)
. (4.23)
Firstly consider counting the number of matrices for a fixed number of elements i, where
r ≤ i ≤ N. Then, there are
|A(i, r)| =
r
∑
q=1
(
i− 1
q− 1
)
(4.24)
matrices having r blocks or fewer, given exacly i elements. From the fact that
i
∑
q=0
(
i
q
)
= 2i, (4.25)
there are 2i−1 matrices in the special case that r = i. To now sum over the number of elements
i, two cases need to be considered, as there cannot be more blocks than elements. The number
of matrices where r ≤ i is given by the summation
N
∑
i=r
|A(i, r)| =
N
∑
i=r
r
∑
q=1
(
i− 1
q− 1
)
. (4.26)
When i < r, there can only be a maximum of i blocks, which means that the number of matrices
is given by
r−1
∑
i=1
|A(i, i)| =
r−1
∑
i=1
2i−1
= 2r−1 − 1. (4.27)
Summing (4.26) and (4.27) gives the required result. 
Theorem 4.6. Given S ∈ A(N, r), |R(S)| = 2N−r.
Proof. As with the proof of Theorem 4.4, represent S as a sequence of block lengths l1, . . . , lr.
Then, by noting that any relaxation of S may be re-blocked to give S, the problem reduces to
finding the number of ways of replacing each ‘’ by either an addition sign (+) or comma (,)
in equation (4.21), subject to r− 1 commas already having been placed (corresponding to S).
This is again a binary decision problem, so the number of ways of placing exactly i commas in
the remaining N − r squares is given by ( N−ri ). The result follows from applying (4.25). 
Lemmas 4.7 - 4.11 are used in the development of a recursively-feasible move-blocked variable
horizon predictive controller.
64
4.3. MOVE BLOCKING FORMULATION
Lemma 4.7. For a blocking matrix S ∈ A(N, r) and some relaxation S′ ∈ R(S, r′), where
N
C
i=1
vi =
(
r
C
q=1
ωq
)
ST, (4.28)
there exists a sequence of block values ω′1, . . . ,ω
′
r′ such that
N
C
i=1
vi =
(
r′
C
q=1
ω′q
)
S′T. (4.29)
Proof. Given S′ ∈ R(S, r′), there exists a blocking matrix M ∈ A(r′, r) such that S′M = S.
Hence
N
C
i=1
vi =
(
r
C
q=1
ωq
)
MTS′T. (4.30)
Then, the result follows by choosing
r′
C
q=1
ω′q =
(
r
C
q=1
ωq
)
MT. (4.31)

Lemma 4.8. For any S ∈ A(N, r) and mapping f : R 7→ Rn×1,
N−1
C
j=0
f (µ(S, j)) =
(
C
q∈S
f (Λ(S, q)− 1)
)
ST. (4.32)
Proof. Firstly, the blocking matrix is decomposed into blocking vectors b1, . . . , br as in (4.9).
Then, it can be seen that(
C
q∈S
f (Λ(S, q)− 1)
)
ST = C
q∈S
f (Λ(S, q)− 1)bTq . (4.33)
As each bTq is a row vector of ones with length `(S, q), the RHS of (4.33) has length N. From the
definition of µ( · , · ) (4.16), expanding the qth term of the RHS gives
f (Λ(S, q)− 1)bTq =
Λ(S,q)−1
C
j=Λ(S,q−1)
f (µ(S, j)). (4.34)
Concatenating over all blocks gives the required result. 
As a numerical example, consider blocking matrix S corresponding to block lengths {3, 2, 2}.
For j = Z[0,6], it is clear that µ(S, j) = {2, 2, 2, 4, 4, 6, 6}. Also, the block indices are given
by q = {1, 2, 3}, so Λ(S, q)− 1 = {2, 4, 6}. It is now easy to see that the result holds when
concatenating the sequences into vectors.
Lemma 4.9. For any S ∈ A(N, r) and S′ ∈ R(S),
µ(S′, j) ≤ µ(S, j), ∀j ∈ Z[0,N−1]. (4.35)
65
4. MOVE BLOCKING FOR COMPLEXITY REDUCTION
Proof. Since any relaxation of a blocking matrix is itself re-blocked to give the original matrix,
consider the qth block of S. In the relaxation, this block is either unchanged or itself divided into
smaller blocks. This means that `(S′, q) cannot be greater than `(S, q), so Λ(S′, q) ≤ Λ(S, q).
The result then follows from (4.16). 
Lemma 4.10. For any S ∈ A(N, r)
Λ(F(S), q) = Λ(S, q + δ[`(S, 1)− 1])− 1, (4.36)
for all q ∈ Z[0,r−δ[`(S,1)−1]].
Proof. Shifting removes the first element of the first block, which reduces the cumulative
block lengths by unity. In addition, when `(S, 1) = 1, the entire first block is removed, so
`(F(S), q) = `(S, q + 1). 
Lemma 4.11. For any blocking matrix S 6= 1,
µ(S, j + 1)− 1 = µ(F(S), j). (4.37)
Proof. The proof follows from (4.16) and Lemma 4.10. It can be seen that
µ(S, j + 1)− 1
=
{
Λ(S, q)− 1− 1,
Λ(S, q− 1)− 1 < j < Λ(S, q)− 1, q ∈N[1+δ[`(S,1)−1],r]
=
{
Λ(F(S), q)− 1,
Λ(F(S), q− 1) < j < Λ(F(S), q), q ∈N[1,r−δ[`(S,1)−1]]
= µ(F(S), j). 
4.4 Blocked VH-MPC
A variable horizon predictive controller will now be formulated using move blocking for
complexity reduction. The variable horizon is achieved by defining a number of fixed horizon
optimisation problems with different horizon lengths and compatible distinct blocking matrices.
Denote the maximum horizon length as N¯ and the maximum allowable number of input
variables in the optimisation problem as r¯. Then, there are |S(N¯, r¯)| possible optimisation
problems corresponding to distinct blocking matrices. A nominal controller will be designed
first using some subset of these, with constraint tightening applied subsequently for robustness
to bounded disturbances.
66
4.4. BLOCKED VH-MPC
4.4.1 Nominal Controller
For a given blocking matrix S ∈ A(NS, rS) ⊆ S(N¯, r¯) (that is, where NS ≤ N¯ and rS ≤ r¯),
define the optimisation problem P(S) as finding
J[P(S)] = min
{υ(S, · )}
NS−1
∑
j=0
(1+ ‖Θ(y(k + j|k)− yr)‖), (4.38)
subject to the dynamics constraints
x(k|k) = x(k) (4.39a)
x(k + j + 1|k) = Ax(k + j|k) + Bu(k + j|k) (4.39b)
y(k + j|k) = Cx(k + j|k) + Du(k + j|k), (4.39c)
the operating constraints
x(k + j|k) ∈ X (4.40a)
u(k + j|k) ∈ U (4.40b)
x(k + NS|k) ∈ T , (4.40c)
and the move blocking constraint
NS−1
C
j=0
u(k + j|k) =
( rS
C
q=1
υ(S, q)
)
ST. (4.41)
The notation υ(S, q) represents the value of the qth block given blocking matrix S and yr is the
output reference defined in Section 4.2.
Remark 4.12. Zero-based indexing is used in the LHS of (4.41) to be consistent with prediction
model notation.
Define now a set of blocking matrices C ⊆ S(N¯, r¯) and a relaxation function
ρ : S(N¯, r¯) 7→ S(N¯, r¯), (4.42)
where
ρ(S) ∈ R(S, r′), r′ ≤ r¯. (4.43)
This function maps a given blocking matrix to precisely one of its relaxations having a max-
imum of r¯ blocks, and is a design choice. Algorithm 4.1 describes the implementation of the
blocked variable horizon in terms of the set C and the function ρ( · ). At each time step, it
essentially solves all problems in C plus the shifted problem, choosing the blocking matrix that
gives the best cost performance.
67
4. MOVE BLOCKING FOR COMPLEXITY REDUCTION
Algorithm 4.1: Nominal Blocked VH-MPC
1 while x(k) /∈ T do
2 if k = 0 then
3 C f ←− {S ∈ C | P(S) feasible};
4 if C f = ∅ then
5 No feasible solution exists for the given starting state and set of blocking matrices;
6 exit;
7 end
8 else
9 Sˆ(k)←− F(S∗′(k)), where S∗′(k) = ρ(S∗(k)) ;
10 C f ←− {S ∈ C | P(S) feasible} ∪ {Sˆ(k)} ;
11 end
12 J∗(x(k))←− min
S∈C f
J[P(S)]; S∗(k)←− arg min
S∈C f
J[P(S)];
13 u(k)←− υ(S∗(k), 1) =⇒ x(k + 1) = Ax(k) + Bυ(S∗(k), 1);
14 k← k + 1;
15 end
Remark 4.13. The set of repeatedly shifted and relaxed blocking matrices
P =
⋃
S∈C
{S, F(ρ(S)), F(ρ(F(ρ(S)))), . . . , 1}, (4.44)
can be calculated in advance if an a priori definition of the problem P(Sˆ( · )) is required for all
time steps.
Remark 4.14. C f is defined as the set of blocking matrices in C whose corresponding optim-
isation problems have a feasible solution, together with the shifted blocking matrix of the
previously optimal problem. It is calculated by actually solving all of the optimisation prob-
lems corresponding to the matrices in C, which implies the solution of multiple optimisation
problems at each time step. However, these problems can be solved in parallel, meaning that
the controller latency is that of the slowest single optimisation problem rather than the sum of
the solution times. Note that either |C| or |C|+ 1 problems are solved at each step, depending
on whether or not it is the first iteration and whether or not the shifted blocking matrix is
already in C.
Theorem 4.15 (Nominal Recursive Feasibility). Given the optimisation problem and associated state
trajectory, input trajectory, output trajectory and horizon length that correspond to the optimal solution
found by Algorithm 4.1 at time k, namely P(S∗(k)), {x∗(k + j|k)}, {u∗(k + j|k)}, {y∗(k + j|k)},
and N∗(k), a feasible solution to the problem P(Sˆ(k)) at time k + 1 is given by
xˆ(k + j + 1|k + 1) = x∗(k + j + 1|k), ∀j ∈ Z[0,N∗(k)−2] (4.45a)
uˆ(k + j + 1|k + 1) = u∗(k + j + 1|k), ∀j ∈ Z[0,N∗(k)−1] (4.45b)
yˆ(k + j + 1|k + 1) = y∗(k + j + 1|k), ∀j ∈ Z[0,N∗(k)−2] (4.45c)
Nˆ(k + 1) = N∗(k)− 1. (4.45d)
68
4.4. BLOCKED VH-MPC
Proof. The candidate solution is evidently feasible with respect to the system dynamics and
constraints, as it is simply a shifted version of the inputs and states at time k, albeit taking
one fewer step to reach completion. To show that the blocking constraints are satisfied by the
candidate inputs, observe that from Lemma 4.7,
N∗(k)−1
C
j=0
u∗(k + j|k) =
(
C
q∈S∗(k)
υ(S∗(k), q)
)
S∗T(k)
=
(
C
q∈S∗′(k)
υ′(S∗′(k), q)
)
S∗′T(k),
(4.46)
where S∗′(k) = ρ(S∗(k)) ∈ A(N∗(k), r′) and υ′( · , · ) are the block values corresponding to the
relaxation. Then,
N∗(k)−2
C
j=0
u∗(k + j + 1|k) =
(
r′
C
q=s+1
υ′(S∗′(k), q)
)
FT(S∗′(k)), (4.47)
where s = δ[`(S∗′(k), 1)− 1]. Hence, substituting the candidate solution gives
Nˆ(k+1)−1
C
j=0
uˆ(k + j + 1|k + 1) =
(
C
q∈Sˆ(k)
υˆ(Sˆ(k), q)
)
SˆT(k), (4.48)
where υˆ(Sˆ(k), q) = υ′(S∗′(k), q + s). This is the required result. 
By induction, therefore, the entire optimisation problem is feasible to completion from an initial
feasible solution. Finite-time completion of the control manœuvre can now be proven.
Theorem 4.16 (Nominal Finite-Time Completion). Algorithm 4.1 guarantees that T will be reached
in finite-time bJ∗(x(0))c, given C f is non-empty at the first iteration.
Proof. It can be seen that the cost of the feasible trajectory at time k + 1, denoted Jˆ( · ) is related
to the cost of the optimal trajectory at time k by the equation
Jˆ(x(k + 1)) = J∗(x(k))− 1. (4.49)
By optimality, this is an upper bound on the optimal cost at time k + 1, therefore
J∗(x(k + 1)) ≤ J∗(x(k))− 1. (4.50)
As the cost function is non-negative by construction, this means that, given an initial feasible
solution and recursive feasibility from Theorem 4.15, completion occurs in at most bJ∗(x(0))c
steps from the starting state. 
69
4. MOVE BLOCKING FOR COMPLEXITY REDUCTION
4.4.2 Robustness
The controller formulation is now extended to include robustness to bounded disturbances.
A disturbance-feedback constraint tightening approach [KRH07] is used to achieve this ro-
bustness. A bounded disturbance is now introduced into the linear system model, giving the
dynamic equation
x(k + 1) = Ax(k) + Bu(k) + w(k), (4.51)
where w(k) ∈ W ⊂ Rn, for some compact setW containing the origin.
Define the tightened input constraints by the recursion
U (0) = U (4.52a)
U (j + 1) = U (j)	 P(j)W . (4.52b)
As in the previous chapter, {P(j) ∈ Rm×n} is a sequence of matrices representing a candidate
control policy providing direct feedback on the disturbance. The policy will be used to ensure
robustness to bounded disturbances, albeit in a modified manner to accommodate move
blocking.
Now, in terms of the function µ( · , · ), define the matrices
L(S, 0) = In (4.53a)
L(S, j + 1) = AL(S, j) + BP(µ(F(ρ(S)), j)) (4.53b)
and the corresponding sets
X (S, 0) = X (4.54a)
T (S, 0) = T (4.54b)
X (S, j + 1) = X (F(ρ(S)), j)	 L(S, j)W (4.54c)
T (S, j + 1) = T (F(ρ(S)), j)	 L(S, j)W . (4.54d)
Remark 4.17. The sets defined by (4.54c) and (4.54d) can be computed recursively by using
the setP (4.44), starting from the blocking matrix S relaxed and shifted to unity, and working
backwards. This tightening scheme is substantially different to that used in [KRH07], as a
different sequence of tightened input, state and terminal constraint sets is needed for each
blocking matrix in the regime.
For a given blocking matrix S ∈ A(NS, rS), define the robust optimisation problem P˜(S) as
finding the minimum (4.38) subject to the tightened constraints
u(k + j|k) ∈ U (µ(S, j)) (4.55a)
x(k + j|k) ∈ X (S, j) (4.55b)
x(k + NS|k) ∈ T (S, NS), (4.55c)
70
4.4. BLOCKED VH-MPC
with the move-blocking constraint (4.41). Note how the function µ( · , · ) is used here to ensure
that the input constraints applied to each block correspond to the tightest constraint for any
input within the block. This is evident from the definition of µ( · , · ) as the block-end index.
By applying Algorithm 4.1 with the nominal problems P( · ) replaced by the robust problems
P˜( · ), recursive feasibility is guaranteed for any allowable sequence of disturbances. This
result is presented in the following theorem.
Theorem 4.18 (Robust Recursive Feasibility). Given the optimisation problem and associated state
trajectory, input trajectory, output trajectory and horizon length that correspond to the optimal solution
found by Algorithm 4.1 at time k, namely P˜(S∗(k)), {x∗(k + j|k)}, {u∗(k + j|k)}, {y∗(k + j|k)} and
N∗(k), a feasible solution to the problem P˜(Sˆ(k)) at time k + 1 is given by
uˆ(k + j + 1|k + 1) = u∗(k + j + 1|k) + P(µ(F(S∗′(k)), j))w(k), ∀j ∈ Z[0,N∗(k)−2] (4.56a)
xˆ(k + j + 1|k + 1) = x∗(k + j + 1|k) + L(S∗(k), j)w(k), ∀j ∈ Z[0,N∗(k)−2] (4.56b)
yˆ(k + j + 1|k + 1) = y∗(k + j + 1|k) + (CL(S∗(k), j) + DP(µ(F(S∗′(k)), j)))w(k),
∀j ∈ Z[0,N∗(k)−2]
(4.56c)
Nˆ(k + 1) = N∗(k)− 1, (4.56d)
where S∗′(k) = ρ(S∗(k)) ∈ A(N∗(k), r′).
Proof. To demonstrate that the candidate solution is feasible, it must be shown to satisfy the
initial constraints, dynamics constraints, move blocking constraints, state and input constraints
and the terminal constraint enforced by problem P˜(Sˆ(k)) = P˜(F(S∗′(k))).
Firstly, to show that the initial constraint holds, observe that the optimal solution satisfies
x∗(k|k) = x(k) (4.57)
from (4.39a). Using (4.39b), (4.41) and substituting the candidate solution (4.56), it can be seen
that
x∗(k + 1|k) = Ax(k) + Bu∗(k|k) (4.58)
=⇒ xˆ(k + 1|k + 1) = Ax(k) + Bυ(S∗(k), 1) + w(k)
= x(k + 1) (4.59)
as required. For the dynamics constraints, the optimal solution satisfies
x∗(k + j + 2|k) = Ax∗(k + j + 1|k) + Bu∗(k + j + 1|k). (4.60)
Substituting the candidate solution gives
xˆ(k + j + 2|k + 1)− L(S∗(k), j + 1)w(k)
= A(xˆ(k + j + 1|k + 1)− L(S∗(k), j)w(k)) + B(uˆ(k + j + 1|k + 1)− P(µ(F(S∗′(k)), j))w(k))
= Axˆ(k + j + 1|k + 1) + Buˆ(k + j + 1|k + 1)− (AL(S∗(k), j) + BP(µ(F(S∗′(k)), j)))w(k)
71
4. MOVE BLOCKING FOR COMPLEXITY REDUCTION
= Axˆ(k + j + 1|k + 1) + Buˆ(k + j + 1|k + 1)− L(S∗(k), j + 1)w(k),
hence,
xˆ(k + j + 2|k + 1) = Axˆ(k + j + 1|k + 1) + Buˆ(k + j + 1|k + 1) (4.61)
as required. For the output dynamics, the optimal solution satisfies
y∗(k + j + 1|k) = Cx∗(k + j + 1|k) + Du∗(k + j + 1|k). (4.62)
Then, substituting the candidate solution shows that
yˆ(k + j + 1|k + 1)− (CL(S∗(k), j) + DP(µ(F(S∗′(k)), j)))w(k)
= C(xˆ(k + j + 1|k + 1)− L(S∗(k), j)w(k)) + D(uˆ(k + j + 1|k + 1)− P(µ(F(S∗′(k)), j))w(k))
= Cxˆ(k + j + 1|k + 1) + Duˆ(k + j + 1|k + 1)− (CL(S∗(k), j) + DP(µ(F(S∗′(k)), j)))w(k).
(4.63)
Hence,
yˆ(k + j + 1|k + 1) = Cxˆ(k + j + 1|k + 1) + Duˆ(k + j + 1|k + 1). (4.64)
For the move blocking constraints, observe that
N∗(k)−2
C
j=0
uˆ(k + j + 1|k + 1) =
N∗(k)−2
C
j=0
u∗(k + j + 1|k) +
N∗(k)−2
C
j=0
P(µ(F(S∗′(k)), j))w(k). (4.65)
Looking at each part of the RHS individually, it is evident that
N∗(k)−2
C
j=0
u∗(k + j + 1|k) =
(
C
q∈Sˆ(k)
υˆ(Sˆ(k), q)
)
SˆT(k), (4.66)
which is shown in the proof of Theorem 4.15. Also, from Lemmas 4.8 and 4.10,
N∗(k)−2
C
j=0
P(µ(F(S∗′(k)), j))w(k) =
(
r′−s
C
q=1
P(Λ(F(S∗′(k)), q)− 1)w(k)
)
FT(S∗′(k)) (4.67)
=
(
C
q∈Sˆ(k)
P(Λ(Sˆ(k), q)− 1)w(k)
)
SˆT(k), (4.68)
where s = δ[`(S∗′(k), 1)− 1]. Combining (4.66) and (4.68) gives
Nˆ(k+1)−1
C
j=0
uˆ(k + j + 1|k + 1) =
(
C
q∈Sˆ(k)
υˆ(Sˆ(k), q) + P(Λ(Sˆ(k), q)− 1)w(k)
)
SˆT(k). (4.69)
Hence, the blocking constraints on the inputs are respected.
For the input constraints, (4.52a)-(4.52b), together with Lemmas 4.9 and 4.11 show that
u∗(k + j + 1|k) ∈ U (µ(S∗(k), j + 1)) (4.70)
72
4.4. BLOCKED VH-MPC
⊆ U (µ(S∗′(k), j + 1)) (4.71)
= U (µ(F(S∗′(k)), j) + 1) (4.72)
= U (µ(F(S∗′(k)), j))	 P(µ(F(S∗′(k)), j))W . (4.73)
This means that, from (2.55),
u∗(k + j + 1|k) + P(µ(F(S∗′(k)), j))w(k) ∈ U (µ(F(S∗′(k)), j))
which implies that
uˆ(k + j + 1|k + 1) ∈ U (µ(F(S∗′(k)), j)) = U (µ(Sˆ(k), j)). (4.74)
For the state constraints
x∗(k + j + 1|k) ∈ X (S∗(k), j + 1). (4.75)
Then from (2.55) and (4.54c),
x∗(k + j + 1|k) + L(S∗, j)w(k) ∈ X (F(S∗′(k)), j) =⇒ xˆ(k + j + 1|k + 1) ∈ X (Sˆ(k), j). (4.76)
Finally, for the terminal constraint
x∗(k + N∗(k)|k) ∈ T (S∗(k), N∗(k)). (4.77)
Then, from (2.55) and (4.54d),
x∗(k + N∗(k)|k) + L(S∗(k), N∗(k)− 1)w(k) ∈ T (F(S∗′(k)), N∗(k)− 1)
=⇒ xˆ(k + 1+ Nˆ(k + 1)|k + 1) ∈ T (Sˆ(k), Nˆ(k + 1)). (4.78)

Theorem 4.19 (Robust Finite-Time Completion). If the cost weightings on the states and inputs are
chosen to ensure that η > 0, where
η = 1− max
w∈W
S∈S(N¯,r¯)
{
‖Θ‖
N¯−2
∑
j=0
∥∥(CL(S, j) + DP(µ(F(S′), j))w∥∥} (4.79)
and S′ = ρ(S), the robust controller guides the system to the target set in at most bJ∗(x(0))/ηc steps,
from the initial state x(0).
Proof. The proof utilises the value function as a Lyapunov-like function. Applying the triangle
inequality and the definition of the induced p-norm, it can be seen that the cost of (4.56),
73
4. MOVE BLOCKING FOR COMPLEXITY REDUCTION
denoted Jˆ(x(k + 1)) is related to the optimal cost at time k, namely J∗(x(k)), by
Jˆ(x(k + 1))− J∗(x(k)) ≤ −1+
N∗(k)−2
∑
j=0
(∥∥Θ(CL(S∗, j) + DP(µ(F(S∗′), j)))w(k)∥∥) (4.80)
≤ max
w∈W
{
N¯−2
∑
j=0
∥∥Θ(CL(S∗, j) + DP(µ(F(S∗′), j)))w∥∥}− 1 (4.81)
≤ max
w∈W
S∈S(N¯,r¯)
{
‖Θ‖
N¯−2
∑
j=0
∥∥(CL(S, j) + DP(µ(F(S′), j)))w∥∥}− 1 (4.82)
= −η. (4.83)
From optimality, Jˆ(x(k + 1)) is an upper bound on the optimal cost at time k + 1, hence
J∗(x(k))− J∗(x(k + 1)) ≥ J∗(x(k))− Jˆ(x(k + 1)) ≥ η. (4.84)
This means that the cost at successive time steps must reduce by at least η. Given that the
cost must be non-negative by construction and η > 0, completion must be reached in at most
bJ∗(x(0))/ηc steps, from initial state x(0). 
4.5 Implementation
The formulation presented in the previous section uses equality constraints on the inputs to
enforce move blocking. This is the most suitable form for analysis, but not necessarily in
implementation. Explicitly implementing the equality constraints in an optimisation program
will not necessarily provide a reduction in computation time. Whether or not a reduction
occurs will depend on how the optimisation algorithm treats equality constraints. In order to
achieve a guaranteed reduction, the prediction model can be condensed to a form that depends
directly on the blocked inputs, removing the equality constraints for move blocking. A further
reduction may be achieved by condensation of the state prediction variables. The condensed
problems defined in the following subsections can be used when implementing Algorithm 4.1.
4.5.1 Input Condensed Prediction Model
Define the input condensed robust optimisation problem P˜C1(S), for S ∈ A(NS, rs) finding the
minimum
J[P˜C1(S)] = min{υ(S, · )}
NS−1
∑
j=0
(1+ ‖Θ(y(k + j|k)− yr)‖), (4.85)
subject to the dynamics constraints
x(k|k) = x(k) (4.86a)
x(k + j + 1|k) = Ax(k + j|k) + Bυ(S, σ(S, j)) (4.86b)
74
4.5. IMPLEMENTATION
y(k + j|k) = Cx(k + j|k) + Dυ(S, σ(S, j)), (4.86c)
and the tightened state and input constraints
υ(S, q) ∈ U (Λ(S, q)− 1) (4.87a)
x(k + j|k) ∈ X (S, j) (4.87b)
x(k + NS|k) ∈ T (S, NS), (4.87c)
where the function σ( · , · ) (4.17) is used to return the block number corresponding to j.
4.5.2 State & Input Condensed Prediction Model
In order to condense the state variables, define parameter-varying matrices by the recursion
A¯(0) = In, A¯(i + 1) = AA¯(i) (4.88)
B¯(0) = 0n×m, B¯(i + 1) = AB¯(i) + B. (4.89)
Define also the condensed state prediction variable for block q by φ(S, q). Then, define the
input and state condensed robust optimisation problem P˜C2(S), for S ∈ A(NS, rs) as finding
the minimum
J[P˜C2(S)] = min{υ(S, · )}
NS−1
∑
j=0
(1+ ‖Θ(y(k + j|k)− yr)‖), (4.90)
subject to the dynamics constraints
φ(S, 1) = x(k) (4.91a)
φ(S, q + 1) = A¯(`(S, q))φ(S, q) + B¯(`(S, q))υ(S, q) (4.91b)
y(k + j|k) = C(A¯(j− t)φ(S, σ(S, j)) + B¯(j− t)υ(S, σ(S, j))) + Dυ(S, σ(S, j)), (4.91c)
and the tightened state and input constraints
υ(S, q) ∈ U (Λ(S, q)− 1) (4.92a)
A¯(i)φ(S, q) + B¯(i)υ(S, q) ∈ X (S,Λ(S, q− 1) + i), ∀i ∈ Z[0,`(S,q)−1] (4.92b)
φ(S, rS + 1) ∈ T (S, NS), (4.92c)
where q ∈N[1,rS] and t = µ(S, j)− `(S, σ(S, j)) + 1 is the block-initial index.
Remark 4.20. It is possible to condense the entire problem, removing all equality constraints.
The resulting problem will then depend only on the initial state and the block values (a so-called
dense MPC formulation). The degree of condensation that gives the greatest computational
benefit will depend on the nature of the system and the solver used, as some solvers can exploit
the sparsity of a problem in order to quicken computation.
75
4. MOVE BLOCKING FOR COMPLEXITY REDUCTION
4.6 Example
From the general controller specification described in the previous sections, there are many
possible ways of choosing the subset of S(N¯, r¯). The simplest method is to subdivide the entire
horizon into equal-length blocks and have the regime consist of matrices corresponding to
multiples of blocks. For some applications, however, it may be beneficial to utilise the flexibility
of the formulation to apply input blocking in only certain parts of the horizon. A particular
method is presented in this section, as may be applicable to vehicle manœuvring problems.
4.6.1 The Partitioned Horizon
The partitioned horizon regime involves blocking matrices where no blocking is applied at
the start and end of the manœuvre, but only in the middle. This is equivalent to partitioning
the prediction horizon into three sections that will be referred to as the initial, medial and
terminal sub-horizons. By applying move blocking only in the medial sub-horizon, a reduction
in complexity is achieved whilst maintaining fine-grained control action at the start and end
of a control manœuvre, as is pertinent to an application like spacecraft rendezvous. This
formulation is inspired by the rubber horizon concept presented in [Har10]. It differs from
the latter by explicitly ensuring pointwise-in-time constraint satisfaction within the medial
sub-horizon.
4.6.2 Formulation
Define Ω1,Ω3 : Ω1 +Ω3 < N¯ as the lengths of the initial and final sub-horizons respectively.
Also define Ω2 = N¯ −Ω1 −Ω3. Assume that a block length of λ is desired in the medial
sub-horizon. The partitioned horizon regime comprises blocking matrices corresponding to
all horizon lengths unblocked up to Ω1 +Ω3 and subsequently adds blocks of length λ to the
medial sub-horizon until the final blocking matrix corresponds to N¯ variables. The final block
added may be shorter than λ if (Ω2 mod λ) 6= 0. The sequence of blocking matrices Si in the
regime can be expressed mathematically as
Si =

Ii 1 ≤ i ≤ Ω1 +Ω3
IΩ1 ⊕ (Ii−Ω1−Ω3 ⊗ 1λ)⊕ IΩ3 Ω1 +Ω3 < i < i¯
IΩ1 ⊕ 1ξ ⊕ (Ii−Ω1−Ω3−1 ⊗ 1λ)⊕ IΩ3 i = i¯,
(4.93)
where i¯ = Ω1 +Ω3 + dΩ2/λe is the total number of blocking matrices in the sequence and
ξ = ((Ω2 − 1) mod λ) + 1 is the length of the first block in the medial sub-horizon when
i = imax. Due to the construction of the regime, i¯ also corresponds to the maximum number of
input variable blocks in any single optimisation problem when the regime is applied.
Then, Algorithm 4.1 is applied to the set of optimisation problems defined by the blocking
matrices Si. The relaxation applied before shifting is defined as follows. When move blocking
76
4.6. EXAMPLE
1 2 3 4 5 6
1 2 4 5 6
1 2 5 6
2 5 6
5 6
6S1
S2
S3
S4
S5
S6
Figure 4.2: Illustration of the partitioned horizon regime
is applied, the blocking matrix S ∈ A(N, r) can be written in the form
S = IΩ1 ⊕ ∆⊕ IΩ3 , (4.94)
where ∆ ∈ A(N −Ω1 −Ω3, r−Ω1 −Ω3) is the blocking matrix corresponding to the move-
blocked part of S. Then, define the relaxation function as
ρ(S) =
IΩ1 ⊕ 1⊕ F(∆)⊕ IΩ3 `(∆, 1) > 1S otherwise. (4.95)
The relaxation unblocks the first element of the medial horizon, which has the effect of only
shifting the medial horizon when F(ρ(S)) is calculated. By using this blocking and relaxation
scheme, Theorem 4.19 ensures that length of the medial sub-horizon in the matrix S∗ reduces
with time, ultimately vanishing. This means that the controller behaves like conventional
VH-MPC formulation after a finite number of steps. Figure 4.2 illustrates the regime for a
simple example, where N¯ = 9, Ω1 = 2 and Ω3 = 2. The structure of ith blocking matrix Si is
shown graphically by the length of the blocks. Note that the medial sub-horizon lies between
the thick lines.
Remark 4.21. As the initial sub-horizon is unblocked, by choosing a constraint-tightening
policy with at most degree-Ω1 nilpotency, the tightened sets corresponding to the states in
the medial and terminal sub-horizons will all be the same. This means that the calculations
described by (4.54c) and (4.54d) simplify significantly, allowing standard constraint tightening
calculations to be used.
4.6.3 Simulation
To illustrate the action of the partitioned horizon, consider a point mass moving in two
dimensions, specified by the 1 Hz discretised dynamics
x(k + 1) = Ax(k) + Bu(k) + w(k), (4.96)
77
4. MOVE BLOCKING FOR COMPLEXITY REDUCTION
where
A =

1 0 1 0
0 1 0 1
0 0 1 0
0 0 0 1
 , B =

0.5 0
0 0.5
1 0
0 1
 , (4.97)
and the disturbance satisfies ‖w(k)‖∞ ≤ 0.45. The control problem is to reach the target set{
x
∣∣∣∣ [200 60 −2 −2]T ≤ x ≤ [202 62 2 2]T} (4.98)
in finite time N, starting from the origin, whilst ideally minimising the cost function
N−1
∑
k=0
1+ 0.3 ‖u(k)‖1 , (4.99)
subject to ‖u(k)‖∞ ≤ 2 and ‖[0 I2]x(k)‖∞ ≤ 5 at all time instants up to completion.
Choosing N¯ = 50 and Ω1 = Ω3 = 5, and using a state and input condensed formulation,
the system is simulated1 for different values of λ over 100 randomly generated admissible
disturbance realisations. A two-step nilpotent policy is used for tightening constraints, which
minimises the number of Pontryagin set difference calculations that are required. The MPC
optimisation problems are solved with a dual-simplex algorithm.
Figure 4.3 shows the effect of different block lengths in the medial sub-horizon on the mean
closed-loop cost. The error bars on the plot indicate one standard deviation. Figure 4.4
illustrates the maximum computation times at each solution time step, for serial and parallel
evaluations of the optimisation problems. The computation times have been normalised to
the maximum serial solution time per time step of the unblocked problem (0.1644 s). Figures
4.5-4.7 show the closed loop trajectories for selected values of λ, whilst Figures 4.8-4.10 show
the input components.
4.6.4 Analysis
It can be seen that the cost is not monotonic with respect to the number of decision variables,
indicating the importance of the blocking structure itself on the closed-loop performance under
disturbances. The fact that the mean cost reduces with reduced complexity of the problem
may at first seem counter-intuitive. However, a possible explanation is the action of the move
blocking, which seems to be affecting the aggressiveness of the control action in responding to
the disturbance. From the control input plots, it is clear that the lower cost is associated with
less aggressive control action being applied. However, after the block lengths are increased
beyond a certain point, the action becomes more aggressive. The maximum latency in both
the serial and parallel cases reduces approximately linearly with a reduction in the number of
1Simulations were performed on an 8-core Mac Pro 3.1 running OS X 10.6.8.
78
4.6. EXAMPLE
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.95
1
1.05
1.1
1.15
N
or
m
al
is
ed
 m
ea
n 
co
st
 to
 c
om
pl
et
io
n
Ratio of maximum number of inputs in blocked to unblocked problem
Figure 4.3: Cost comparison of different block lengths in the partitioned horizon
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Ratio of maximum number of inputs in blocked to unblocked problem
N
or
m
al
is
ed
 m
ax
 s
ol
ut
io
n 
tim
e 
pe
r t
im
es
te
p
 
 
Serial
Parallel
Figure 4.4: Compuation time comparison different block lengths in the partitioned horizon
79
4. MOVE BLOCKING FOR COMPLEXITY REDUCTION
0 20 40 60 80 100 120 140 160 180 200
0
10
20
30
40
50
60
x
y
Figure 4.5: Closed-loop trajectories for λ = 1
0 20 40 60 80 100 120 140 160 180 200
0
10
20
30
40
50
60
x
y
Figure 4.6: Closed-loop trajectories for λ = 4
0 20 40 60 80 100 120 140 160 180 200
0
10
20
30
40
50
60
x
y
Figure 4.7: Closed-loop trajectories for λ = 40
80
4.6. EXAMPLE
0 5 10 15 20 25 30 35
−3
−2
−1
0
1
2
3
4
k
in
pu
t
 
 
input x
input y
Figure 4.8: Control inputs for λ = 1
0 5 10 15 20 25 30 35
−3
−2
−1
0
1
2
3
4
k
in
pu
t
 
 
input x
input y
Figure 4.9: Control inputs for λ = 4
0 5 10 15 20 25 30 35 40 45
−3
−2
−1
0
1
2
3
4
k
in
pu
t
 
 
input x
input y
Figure 4.10: Control inputs for λ = 40
81
4. MOVE BLOCKING FOR COMPLEXITY REDUCTION
input variables. These latencies give an indication of what the maximum sampling rate of the
system can be when considering real-time implementation.
4.7 Conclusion
This chapter has developed a new robust formulation of VH-MPC using move blocking to
reduce computational complexity. The formulation has been shown to preserve the guarantees
of recursive feasibility and finite-time completion in the presence of bounded disturbances.
Using a particular choice of blocking regime, namely the partitioned horizon, the formulation
has been demonstrated on an example system.
4.8 Acknowledgements
The work in this chapter has been published in the journal Systems & Control Letters [SM12].
The simulations in this chapter were carried out using Matlab R2009b, YALMIP [Löf04], GLPK
[Glp] and the Invariant Set Toolbox [Ker00].
82
Chapter 5
Optimal Constraint Tightening Policies
5.1 Introduction
As constraint tightening achieves robustness by restricting the size of constraint sets at future
steps into the prediction horizon, this will result in a reduction in the size of the controller’s
Region of Attraction (ROA), which is the set of states for which an initial feasible solution exists
to the MPC optimisation problem. In the formulations of constraint tightening presented in
[RH03], [RH06a], no restriction is placed on the policy for the guarantees of robust recursive
feasibility and finite-time completion to hold in the closed loop. Of course, it is necessary that
the chosen policy will not cause any of the tightened sets to become empty, which would result
in the controller’s ROA being empty. For convenience in calculating tightened sets, nilpotent
policies have been proposed; for this reason, such policies have been used in the preceding
chapters. However, these policies may not be the best choice if a large size ROA is desired; the
degree of freedom in choosing the CT policy suggests that optimisation is possible.
Kuwata, Richards and How [KRH07] demonstrate that a disturbance feedback parameterisation
allows the constraint tightening policy to form part of a convex optimisation problem. The
problem aims to find the constraint tightening policy able to tolerate strong disturbances, for
a fixed prediction horizon. Essentially it maximises a scaling factor on the disturbance set,
whilst still ensuring non-emptiness of the tightened constraint sets, as well as the existence
of a non-empty robust control-invariant terminal set. A numerical example shows how the
resulting policy attenuates the size of the region of attraction less severely than policies selected
in other ways, as the disturbance set is dilated.
This chapter builds upon the work of [KRH07], but considers a fundamentally different
criterion. It demonstrates how, given a fixed polytopic disturbance set, constraint-tightening
policies parameterised in terms of disturbance feedback can be optimised directly to increase the
volume of the controller’s ROA. Furthermore, this chapter analyses the variable horizon case, so
there is no requirement on terminal set invariance. The variable horizon does, however, require
approximations to be made, given the non-convexity of the problem. An inner ellipsoidal
approximation to a fixed horizon ROA is maximised, which is an inner approximation to the
83
5. OPTIMAL CONSTRAINT TIGHTENING POLICIES
non-convex (in general) union of convex sets that defines the full ROA over all horizon lengths
(denoted the Variable Horizon ROA or VH-ROA). This inner approximation is shown to be
effective in the numerical examples that are used to demonstrate the technique.
The examples show that, compared to using nilpotent policies or other heuristics, the policy
resulting from the optimisation method presented in this chapter produces the largest domain
of attraction volume. Furthermore, it is shown that this policy can drastically reduce the
horizon length required for the ROA to contain a given subset of the state space. The effect of
the policy on closed-loop cost performance is also investigated.
5.1.1 Nomenclature
For any matrix A (including vectors), the notation ATi represents its ith row. The hypervolume
of a given polytope C, or its Lebesgue measure, is denoted vol(C). The projection of a polyhed-
ron A ⊂ Rn onto the subspace Rm, where m < n is denoted by
Pm(A) = {x ∈ Rm | ∃u ∈ R(n−m) : (x, u) ∈ A}. (5.1)
A p-norm ball in n dimensions with radius e is defined by the set
Bnp(e) = {x ∈ Rn | ‖x‖p ≤ e}. (5.2)
5.2 Problem Formulation
As in the last two chapters, consider a discrete-time linear system with states x(k) ∈ Rn, inputs
u(k) ∈ Rm and outputs y(k) ∈ Rp specified by the difference equation
x(k + 1) = Ax(k) + Bu(k) + w(k) (5.3)
y(k) = Cx(k) + Du(k), (5.4)
where A ∈ Rn×n, B ∈ Rn×m, C ∈ Rp×n and D ∈ Rp×m. It is assumed that the pair (A, B) is
controllable. The system is subject to pointwise-in-time output constraints of the form
y(k) ∈ Y = {y | Ey ≤ f , E ∈ Rq×p, f ∈ Rq}, (5.5)
where Y is a polytope defined by q individual half-space constraints. In terms of this set, define
the output admissible set
O = Pn ({(x, u) : E(Cx + Du) ≤ f }) . (5.6)
The control objective is to steer the state to the target set
T = {x | Gx ≤ h, G ∈ Rr×n, h ∈ Rr}, (5.7)
84
5.2. PROBLEM FORMULATION
which is a polytope comprising r individual half-space constraints. This set is to be reached in
N < ∞ steps, expressed as the constraint
x(N) ∈ T . (5.8)
The polytopically-bounded state disturbance satisfies
w(k) ∈ W = {w | ζw ≤ θ}, (5.9)
for a given matrix and vector ζ ∈ Ra×n and θ ∈ Ra respectively. This set is assumed to contain
the origin. The control problem is to steer the system state to T in finite time, whilst minimising
a given cost function.
Remark 5.1. In this chapter, the output constraints are assumed to be polytopic to ensure
that the ROA is compact. If this assumption were not made, certain system dynamics (e.g.
open-loop nilpotent systems) could result in an unbounded ROA.
A disturbance-feedback constraint-tightened predictive controller is designed to address the
control problem. Recursively define the matrices
L(0) = In (5.10a)
L(j + 1) = AL(j) + BP(j) (5.10b)
Q(j) = CL(j) + DP(j) (5.10c)
and the corresponding sets
Y(0) = Y (5.11a)
T (0) = T (5.11b)
Y(j + 1) = Y(j)	Q(j)W (5.11c)
T (j + 1) = T (j)	 L(j)W . (5.11d)
As in previous chapters, the sequence of matrices
P = {P(j) ∈ Rn×m}, (5.12)
is the constraint tightening policy providing direct feedback on the disturbance. The prediction
model is then defined by the dynamics
x(k + j + 1|k) = Ax(k + j|k) + Bu(k + j|k) (5.13a)
y(k + j|k) = Cx(k + j|k) + Du(k + j|k), (5.13b)
subject to the constraints
y(k + j|k) ∈ Y(j) (5.14a)
85
5. OPTIMAL CONSTRAINT TIGHTENING POLICIES
x(k + N(k)|k) ∈ T (N(k)). (5.14b)
5.3 Affine Parameterisation of Tightened Sets
In order to optimise over the disturbance feedback matrices, the Pontryagin difference opera-
tions (5.11c)-(5.11d) need to be parameterised as affine functions of the policy matrices P( · ). In
order to accomplish this, define the sequences of vectors {s(j) ∈ Rq}N−1j=0 and {σ(j) ∈ Rr}Nj=0,
as well as the sequences of matrix variables {Z(j) ∈ Ra×q}N−1j=0 and {Z¯(j) ∈ Ra×r}Nj=0, for some
fixed horizon length N.
Theorem 5.2. For any output y and state x, for all j ∈ Z[0,N−1],
y ∈ Y(j)⇐⇒ ∃{Z(l) ≥ 0}j−1l=0 : Ey ≤ f − s(j), (5.15)
and for all j ∈ Z[0,N]
x ∈ T (j)⇐⇒ ∃{Z¯(l) ≥ 0}j−1l=0 : Gx ≤ h− σ(j), (5.16)
where
s(0) = 0q×1 (5.17a)
s(j + 1) = s(j) + ZT(j)θ (5.17b)
EQ(j) = ZT(j)ζ (5.17c)
and
σ(0) = 0r×1 (5.18a)
σ(j + 1) = σ(j) + Z¯T(j)θ (5.18b)
GL(j) = Z¯T(j)ζ. (5.18c)
Proof. The theorem is proved by induction on j for the tightened output sets. The proof for the
terminal sets is completely analogous. Then, for the base case, it is clear that
y ∈ Y(0)⇐⇒ Ey ≤ f − s(0) = f , (5.19)
as s(0) = 0q×1 and there is no dependence on Z( · ). Now, assume that for any j,
y ∈ Y(j)⇐⇒ ∃{Z(l) ≥ 0}j−1l=0 : Ey ≤ f − s(j), (5.20)
where the constraints (5.17) are satisfied. Then, using (5.11c) and the definition of the Pontry-
agin difference (2.55),
y ∈ Y(j + 1)⇐⇒ y + b ∈ Y(j), ∀b ∈ Q(j)W (5.21)
86
5.3. AFFINE PARAMETERISATION OF TIGHTENED SETS
which, using the induction assumption, means that there exists Z(0), . . . , Z(j− 1) ≥ 0 satisfying
(5.17) such that
E(y + b) ≤ f − s(j), ∀b ∈ Q(j)W (5.22)
Ey + EQ(j)w ≤ f − s(j), ∀w ∈ W . (5.23)
The universal quantifier can then be eliminated to give
Ey +max
w∈W
EQ(j)w ≤ f − s(j) (5.24)
where the maximisation is taken row-wise. Given that each row-wise maximisation defines a
linear program, its (asymmetric) dual problem can be exploited as in [GKM06]. Hence, for the
ith row,
max
w∈W
ETi Q(j)w = minzi
zTi θ (5.25)
subject to
ETi Q(j) = z
T
i ζ (5.26)
zi ≥ 0, (5.27)
for dual variable zi ∈ R1×a. The minimisation can then be removed by observing that, for the
ith row of (5.24),
ETi y +minzi
zTi θ ≤ fi − si(j)⇐⇒ ∃zi ≥ 0 : ETi y + zTi θ ≤ fi − si(j). (5.28)
Applying the same reasoning to all i ∈N[1,q], the dual variables can be combined to give
Ey +max
w∈W
EQ(j)w ≤ f − s(j)⇐⇒ ∃Z(j) ≥ 0 : Ey + ZT(j)θ ≤ f − s(j), (5.29)
where
Z(j) =
q
C
i=1
zi (5.30)
EQ(j) = ZT(j)ζ. (5.31)
Hence, there exists Z(0), . . . , Z(j) ≥ 0 satisfying (5.17) such that
Ey ≤ f − (s(j) + ZT(j)θ) (5.32)
≤ f − s(j + 1). (5.33)
which demonstrates that
y ∈ Y(j + 1)⇐⇒ ∃{Z(l)}jl=0 : Ey ≤ f − s(j + 1), (5.34)
87
5. OPTIMAL CONSTRAINT TIGHTENING POLICIES
if the conditions (5.17) are satisfied. The same method can be used to show the result for the
tightened terminal set constraints, so the proof is omitted for brevity. 
This parameterisation goes further than [KRH07], where the tightened output set Y(N − 1) is
expressed in terms of a series of affine inequalities based on the vertices of the disturbance set,
which can become computationally burdensome for higher-dimensional polytopes. Further-
more, the assumption 0 ∈ Y is made in [KRH07], whereas no such assumption is required in
Theorem 5.2.
Corollary 5.3 shows that the RHS of the polytopes defining the output and terminal constraints
in Theorem 5.2 can never increase, which means that the constraint sets can never grow. This is
expected, given the action of the repeated P-difference operations. The result also verifies the
fact that the sets remain compact under tightening.
Corollary 5.3. In the affine parameterisation presented in Theorem 5.2, the variables s( · ) and σ( · )
are always monotonically increasing in j.
Proof. This result arises from the fact that the disturbance set must contain the origin and the
matrices Z( · ) and Z¯( · ) are element-wise non-negative. The former condition implies that
θ ≥ 0, which, combined with the latter condition ensures that ZT(j)θ ≥ 0 and Z¯T(j)θ ≥ 0, for
all j. These conditions can then be combined with recurrence relationships for s( · ) and σ( · )
defined in (5.17) - (5.18) to show that s( · ) and σ( · ) increase monotonically. 
5.4 Maximising the Region of Attraction
Having defined the tightened sets as affine functions of the constraint tightening policy, the
regions of attraction for fixed horizon lengths can also be parameterised in terms of this policy.
This then allows a measure of the volume of this region to be optimised.
5.4.1 Parameterisation
Defining the block matrix
Φ =

EC ED 0 · · · 0 0
ECA ECB ED · · · 0 0
ECA2 ECAB ECB · · · 0 0
...
...
...
. . .
...
...
ECAN−1 ECAN−2B ECAN−3B · · · ECB ED
GAN GAN−1B GAN−2B · · · GAB GB

, (5.35)
88
5.4. MAXIMISING THE REGION OF ATTRACTION
and vectors
v =
[
x(0)
u
]
, u =

u(0|0)
u(1|0)
...
u(N − 1|0)
 ,ψ =

f
f
...
f
h

, s =

s(0)
s(1)
...
s(N − 1)
σ(N)

, (5.36)
the following lemma shows how the tightened operating constraints can be condensed into
half-space constraints on the initial states and sequence of inputs.
Lemma 5.4. For a given fixed horizon length N, the constraints (5.14) on the initial MPC optimisation
problem can be expressed in the form
Φv ≤ ψ− s, (5.37)
∃{Z(j) ≥ 0a×q}, {Z¯(j) ≥ 0a×r} :
s satisfies (5.17) and (5.18)
(5.38)
Proof. This result is obtained by condensing the states in the dynamics equations (5.13) at time
step k = 0 and enforcing the affine constraints derived in Theorem 5.2 on the resulting outputs.
Specifically,
y(0|0) ∈ Y(0)⇐⇒ E(Cx(0) + Du(0|0)) ≤ f − s(0)
y(1|0) ∈ Y(1)⇐⇒ E(C(Ax(0) + Bu(0|0)) + Du(1|0) ≤ f − s(1)
y(2|0) ∈ Y(2)⇐⇒ E(C(A2x(0) + ABu(0|0) + Bu(1|0)) + Du(2|0)) ≤ f − s(2)
...
y(N − 1|0) ∈ Y(N − 1)⇐⇒ E(C(AN−1x(0) + AN−2Bu(0|0) + · · ·+ Bu(N − 2|0))
+ Du(N − 1|0) ≤ f − s(N − 1)
x(N|0) ∈ T (N)⇐⇒ G(ANx(0) + AN−1Bu(0|0) + · · ·+ Bu(N − 1|0)) ≤ h− σ(N),
where s satisfies (5.17) and (5.18). 
In terms of this constraint set, Theorem 5.5 quantifies the fixed horizon region of attraction for
a given horizon length N.
Theorem 5.5. The set of states R(N,P) for which an initial N-step feasible solution exists, for CT
policy P is given by
R(N,P) =
x
∣∣∣∣∣∣∣∣
Φˆx ≤ ψˆ− Γs
∃{Z(j) ≥ 0a×q}, {Z¯(j) ≥ 0a×r} :
s satisfies (5.17) and (5.18)
 , (5.39)
where the matrices Γ ∈ Rt×(d−n) and Φˆ ∈ Rt×n, as well as vector ψˆ ∈ R define the projection of the
89
5. OPTIMAL CONSTRAINT TIGHTENING POLICIES
polyhedron
A =
{
(s, v)
∣∣∣∣∣ [I Φ]
[
s
v
]
≤ ψ
}
(5.40)
onto the subspace Rd, where
d = (N − 1)q + r + n. (5.41)
That is, {
(s, x(0)) ∈ Rd
∣∣∣∣∣ [Γ Φˆ]
[
s
x(0)
]
≤ ψˆ
}
, Pd(A). (5.42)
Proof. The constraint (5.37) derived in Lemma 5.4 can be rewritten in the form
Φ
[
x(0)
u
]
+ s ≤ ψ (5.43)
=⇒
[
I Φ
]  s[x(0)
u
] ≤ ψ (5.44)
which defines the polyhedron A. By taking its projection onto Rd, what remains is a poly-
hedron defined by inequalities depending only on the initial state x(0) and the vector s. This
polyhedron can be written in the form
[
Γ Φˆ
] [ s
x(0)
]
≤ ψˆ (5.45)
=⇒ Φˆx(0) ≤ ψˆ− Γs. (5.46)
Combining these constraints with the requirements on the dual variables gives the desired
result. 
Remark 5.6. It is assumed that |P| ≥ N − 1. Any extra disturbance feedback gains are simply
unused in the calculation of s( · ) and σ( · ).
The projection operation in Theorem 5.5 can be implemented with the conceptually simple
Fourier-Motzkin elimination technique [DE73], or alternatively, a less computationally intens-
ive, but conceptually complex algorithm such as Equality Set Projection (ESP) [JKM04]. Using
ESP will require an artificial bounding of the polyhedron A to convert it into a polytope. This
can be achieved by adding extra constraints of the form s < M1(d−u), as these are the only
unbounded variables in the absence of the dual constraints. M needs to be chosen large enough
that it does not interfere with the other constraints on the dual variables. A sufficient condition
is to choose M ≥ max{diam(Y), diam(T )}. Regardless of the projection method used, the
computational burden is only incurred offline. Note that it is possible to project away the dual
variables as well, however, in the interests of reducing computational complexity, they are
retained to be defined explicitly in subsequent optimisation problems.
90
5.4. MAXIMISING THE REGION OF ATTRACTION
5.4.2 Optimisation
The volume of a general polytope has no closed-form expression, so instead, optimisation
focuses on maximising the volume of an inner ellipsoidal approximation toR(N,P). Consider
an ellipsoid specified by
E(c,Ω) = {x ∈ Rn | (x− c)T(ΩTΩ)−1(x− c) ≤ 1}, (5.47)
with centre c ∈ Rn and shape matrix Ω ∈ Rn×n. A common assumption is that Ω  0 [ZG03].
Then, Theorem 5.7 states the conditions under which this ellipsoid is contained within the
region of attraction.
Theorem 5.7. E(c,Ω) ⊂ R(N,P) iff
Φˆc + γ(E) ≤ ψˆ− Γs (5.48a)
∃{Z(j) ≥ 0a×q}, {Z¯(j) ≥ 0a×r} :
s satisfies (5.17) and (5.18)
(5.48b)
where
γ(E) =
(
t
C
i=1
‖ΩΦˆi‖2
)T
∈ Rt. (5.49)
Proof. The proof uses the same method as in [ZG03], utilising the parameterised polytope
defining the N-step ROAR(N,P). Firstly, E( · , · ) is re-expressed as an affine transformation
on a unit hypersphere centred at the origin, giving the expression
E(c,Ω) = {x ∈ Rn | x = Ωξ + c : ξ ∈ Rn, ‖ξ‖2 ≤ 1}. (5.50)
Then, it can be seen that all points within the ellipsoid satisfy the half-space constraint defined
by the ith row of (5.46) iff
ΦˆTi (c +Ωξ) ≤ ψˆ− Γs, ∀ξ : ‖ξ‖2 = 1 (5.51)
As in (5.24), the universal quantifier can be replaced by
ΦˆTi c + sup
‖ξ‖2=1
ΦˆTi Ωξ ≤ ψˆi − Γis. (5.52)
Given that the magnitude of ξ is fixed, the supremum is achieved when vectors ΩΦˆi and ξ are
parallel. Hence
ΦˆTi c + ‖ΩΦˆi‖2 ≤ ψˆi − Γis. (5.53)
Defining γ( · ) as in (5.49), the inequalities for all rows can be expressed in the form (5.48a).
The result then follows by enforcing the conditions on the dual variables. 
Theorem 5.8. For a given N, the volume of an inner ellipsoidal approximation to R(N,P) can be
91
5. OPTIMAL CONSTRAINT TIGHTENING POLICIES
maximised by solving the convex constrained optimisation problem
P(N) =
 minc,Ω,P − log detΩsubject to (5.48) and Ω  0 (5.54)
Proof. Given that the volume of an ellipsoid is proportional to the product of the lengths of
its semiaxes, which correspond to the eigenvalues of Ω, the volume is proportional to det(Ω).
The problem is made convex by taking logarithms; it is easy to show that the volume of the
ellipsoid increases monotonically with log det(Ω). Constrained maximisation of this objective
is a Semidefinite Program (SDP), specifically, a convex generalised eigenvalue problem, where
all of the constraints (5.48) can be expressed in LMI form [Boy+94; BV04]. 
Remark 5.9. If it is desirable to maximise the volume of the ROA around some given region
G ⊆ O, then the constraint
c ∈ G (5.55)
can be enforced. In addition, if a nilpotent policy is desired for computational reasons, then
constraints of the form
L(η) = 0n×n
P(η + j) = 0m×n, ∀j ≥ 0
(5.56)
can be enforced, for degree-η nilpotence. These constraints also reduce the complexity of the
projection, as s(η + j) = s(η) and σ(η + j) = σ(η), for all j > 0, allowing the dimension of
s to be reduced. This will require the identity matrix in (5.44) to be replaced appropriately.
It should be noted that enforcing either (5.55) or (5.56) can make it harder to find a feasible
solution to the problem.
The constraint tightening policy that minimises (5.54), denoted P∗(N) defines the policy that
maximises an approximation to the volume of the N-step region of attraction. To use this in a
variable horizon framework, where the maximum horizon length is N¯, it remains to solve all
problemsP(1), . . . ,P(N¯) and choose the policy with horizon length given by
N∗ = arg min
N∈N[1,N¯]
P(N) (5.57)
where, with an abuse of notation, minP( · ) denotes the minimum value of the problem’s
objective function. Then, the policy Pmax-vol = P∗(N∗) will give the maximum ellipsoidal
volume. Once the optimal policy has been found, the tightened sets do not need to be re-
calculated by repeated application of the P-difference operation. This is due to the fact that
the values of s( · ) and σ( · ) that correspond to the optimal policy are already known after
solving (5.54) and hence, can be subtracted from f and h respectively to immediately give the
half-space representation of the tightened sets.
92
5.5. EXAMPLES
For any horizon length, the resulting policy ensures non-emptiness of the tightened sets, as
shown in Corollary 5.10.
Corollary 5.10. The optimal policy P∗(N), if it exists, ensures that the tightened sets
{Y(j)}N−1j=0 , (5.58)
as well as the terminal set T (N), are non-empty.
Proof. This can be easily seen by observing that for every x(0) ∈ R(N,P∗(N)), there exists a
sequence of inputs u such that, for any j ∈N[1,N−1],
y(j|0) = Cx(j|0) + Du(j|0) ∈ Y(j). (5.59)
This input sequence also guarantees the existence of x(N|0) that satisfies
x(N|0) ∈ T (N). (5.60)
Hence, these tightened sets must be non-empty by construction. 
Remark 5.11. Corollary 5.10 shows that, unlike in [KRH07], additional constraints do not need
to be applied to ensure the non-emptiness of the state constraint sets and the terminal set
when the policy is chosen according to (5.54). This formulation does not, however, guarantee
that the sets are of nonzero measure. If such a guarantee is required for improving numerical
conditioning, a sufficient condition is to ensure that there exist ex, ey > 0 such that B
p
p(ex) ⊂
Y(N − 1) and Bnp(ey) ⊂ T (N), for some p-norm balls. Specifically, the region of attraction
(5.39) is restricted with the additional constraints
∃y, x : (5.61a)
Ey + ey
( q
C
i=1
‖Ei‖
)T
≤ f − s(N − 1)
Gx + ex
( q
C
i=1
‖Gi‖
)T
≤ h− σ(N).
(5.61b)
The derivation of these constraints is similar to (5.53).
5.5 Examples
This section presents example scenarios that demonstrate the efficacy of designing a constraint-
tightening policy in the manner detailed in the previous sections. A two-state system is
considered first, which allows the region of attraction to be easily visualised and its exact
volume to be readily computed. This is used to show the superiority of the max-vol policy by
comparing actual ROA volumes over different policies. In addition, the effect of the policy on
93
5. OPTIMAL CONSTRAINT TIGHTENING POLICIES
closed-loop costs is also analysed. A final scenario presents a three-state system which further
demonstrates the benefit of using max-vol as opposed to using a simple heuristic.
5.5.1 Symmetric Target Region Placement
Consider the unstable discrete-time linear system defined by state-space matrices
A =
[
1.2 1.2
0 1
]
, B =
[
0.5
1
]
, C =
1 00 1
0 0
 , D =
00
1
 (5.62)
The outputs of the system comprise the states and the inputs fed through, that is, y = (x, u).
Given a maximum prediction horizon length of N¯ = 4 steps, the control objective is to steer
the system to the target set
T1 =
{
x
∣∣∣∣ [−2 −2]T ≤ x ≤ [2 2]T} , (5.63)
whilst satisfying the output constraints
Y =
{
y
∣∣∣∣ [−15 −5 −3]T ≤ y ≤ [15 5 3]T} (5.64)
and being robust to the bounded disturbance
W =
{
w
∣∣∣∣ [−0.3 −1]T ≤ w ≤ [0.3 1]T} . (5.65)
This is a problem with symmetric constraints and a target set that is placed centrally within
the output-admissible set. Solving (5.57) shows that the maximum ellipsoid volume results
from choosing N = 4, corresponding to the maximum horizon length. Figure 5.1 illustrates
the nominal (untightened) region of attraction as well as the 4-step ROA corresponding to
Pmax-vol. The terminal set and inner ellipsoidal approximation are also shown. Figure 5.2
shows the regions of attraction for all horizon lengths using the same policy. Note that the
1-step ROA is policy-independent, since CT over a prediction horizon of one tightens the
terminal set directly by W , with no input tightening applied. It is evident that the ROA
corresponding to the maximum horizon length contains all of the others (apart from N = 1
which is policy-independent), suggesting that the approximation is reasonable for this system.
Table 5.1 compares the exact volumes for the 4-step ROA , as well as the total volume over all
horizon lengths, using different tightening policies. These ROAs are illustrated in Figure 5.3.
Exact volume calculations are achieved using Delaunay triangulation [Del34] for the convex
sets and in the case of the non-convex union over horizon lengths, the inclusion-exclusion
principle is applied [Com74]. The volumes for the entry marked “nominal” correspond to those
for an arbitrary policy whenW = ∅ i.e. the untightened non-robust problem. In addition, the
94
5.5. EXAMPLES
policy labelled “max-dist” is an implementation of the policy designed in [KRH07], adapted to
the variable horizon problem. The derivation of this policy can be found in Appendix A. Plots
of the full VH-ROAs for all policies can also be found there.
These results demonstrate the superiority of the max-vol policy for this application, as far
as the size of the ROA is concerned. The policy provides a 36.6% larger 4-step ROA and a
35.6% larger VH-ROA as compared to max-dist policy. This is not surprising however, as
the max-dist strategy poses a fundamentally different optimisation problem, as discussed
in Section 5.1. The nilpotent policies are seen to drastically reduce the region of attraction,
suggesting the importance of optimising the constraint tightening policy if a large ROA is
desired. Subsection 5.5.2 vividly illustrates this point.
P vol (R(4,P)) vol
(⋃4
i=1R(i,P)
)
nominal 199.2115 199.2115
max-vol 93.8338 94.0955
max-dist 59.4968 60.6070
3-step nilpotent 39.4018 40.8979
2-step nilpotent 23.2753 28.0190
Table 5.1: ROA volumes over different policies for T1
x1
x 2
-15 -10 -5 0 5 10 15
-5
-4
-3
-2
-1
0
1
2
3
4
5
Nominal
Max-vol
Terminal set
Inner ellipsoid
Figure 5.1: Four-step ROA for T1 using Pmax-vol
95
5. OPTIMAL CONSTRAINT TIGHTENING POLICIES
x1
x 2
-8 -6 -4 -2 0 2 4 6 8
-5
-4
-3
-2
-1
0
1
2
3
4
5
N=4
N=3
N=2
N=1
Terminal set
Figure 5.2: Total variable horizon ROA for T1 using Pmax-vol
x1
x 2
-15 -10 -5 0 5 10 15
-5
-4
-3
-2
-1
0
1
2
3
4
5
Nominal
Max-vol
Max-dist
3-step Nil
2-step Nil
Terminal set
Figure 5.3: Four-step ROA comparison over different policies for T1
96
5.5. EXAMPLES
5.5.2 Horizon Length Reduction
Choosing the max-vol policy can provide a drastic reduction in the required maximum horizon
length, directly decreasing online computational complexity. Considering once again the
control problem posed in Subsection 5.5.1, Figure 5.4 shows the 4-step ROA for Pmax-vol
compared to the 14-step ROA for a two-step nilpotent tightening policy. For both policies,
the full variable horizon ROA, apart from the 1-step ROA which is policy-independent, is
contained within the 4-step ROA. It is evident that even with a horizon length over three
times as long, the nilpotent policy is unable to give the region of attraction volume achievable
by using Pmax-vol. This suggests that simply choosing a 2-step nilpotent policy for ease of
computing the P-difference as in [RH06a] can significantly increase the required horizon length
if the ROA is required to contain a given subset of O that is far from the target set.
x1
x 2
-8 -6 -4 -2 0 2 4 6 8
-5
-4
-3
-2
-1
0
1
2
3
4
5
4-step ROA for max-vol
14-step ROA for 2-step nil
Terminal set
Figure 5.4: Comparison of 4-step ROA for T1 using Pmax-vol with 14-step ROA for T1 using
P2-step-nil
5.5.3 Asymmetric Target Region Placement
The effect of allowing a variable centre for the ellipsoid in (5.54) can be illustrated by shifting
the terminal set. Using the same system dynamics, the new target region is given by
T2 =
{
x
∣∣∣∣ [11 1]T ≤ x ≤ [15 5]T} , (5.66)
97
5. OPTIMAL CONSTRAINT TIGHTENING POLICIES
which is the upper-right corner of the output-admissible set. It can be seen in Figure 5.5 that
the ellipsoid centre has moved to accommodate the asymmetry. The effect of this new target is
especially evident when looking at the total variable horizon ROA for the max-vol policy, as
illustrated in Figure 5.6. The 4-step ROA no longer contains all of the other ROAs for smaller
horizon lengths. Table 5.2 compares the volumes as before, where it is seen that max-vol still
provides the largest volume of all of the policies. However, the increase in VH-ROA volume
over the max-dist policy is more modest (14%). As in the previous subsection, the 4-step ROAs
for all policies are illustrated in Figure 5.7.
P vol (R(4,P)) vol
(⋃4
i=1R(i,P)
)
nominal 167.0490 174.0646
max-vol 58.5689 99.0855
max-dist 34.4009 86.7941
3-step nilpotent 24.5657 80.5391
2-step nilpotent 17.3852 74.1435
Table 5.2: ROA volumes over different policies for T2
x1
x 2
-10 -5 0 5 10
-5
-4
-3
-2
-1
0
1
2
3
4
5
Nominal
Max-vol
Terminal set
Inner ellipsoid
Figure 5.5: Four-step ROA for T2 using Pmax-vol compared with nominal ROA
98
5.5. EXAMPLES
x1
x 2
-5 0 5 10 15
-4
-3
-2
-1
0
1
2
3
4
5
N=4
N=3
N=2
N=1
Terminal set
Figure 5.6: Total variable horizon ROA for T2 using Pmax-vol
x1
x 2
-10 -5 0 5 10
-5
-4
-3
-2
-1
0
1
2
3
4
5
Nominal
Max-vol
Max-dist
3-step Nil
2-step Nil
Terminal set
Figure 5.7: Four-step ROA comparison over different policies for T2
99
5. OPTIMAL CONSTRAINT TIGHTENING POLICIES
5.5.4 Effect on Cost
It is possible to evaluate the cost performance of different tightening policies by choosing initial
states randomly from a starting region that is contained within the 4-step regions of attraction
of all of the policies considered. Considering once again the problem of steering the system
state to T2, a suitable starting region is given by
X0 =
{
x
∣∣∣∣ [−3 2]T ≤ x ≤ [−2 3]T} . (5.67)
The cost function
N−1
∑
j=0
1+ ‖Θ(y(k + j|k)− yr)‖1 , (5.68)
where
Θ = diag
{[
0.1 0.1 1
]}
, yr =
[
0 13 3
]
, (5.69)
penalises a weighted sum of time to completion, input effort, and distance from the centre
of T2. Table 5.3 shows the mean cost over 100 random starting states within X0 and random
admissible disturbance sequences, for different tightening policies. It is evident that optimising
the ROA volume does not directly correspond to cost optimality, as the max-dist policy gives
the best cost performance. It should be noted, however, that the differences in cost are minimal
in comparison to the differences in volume of the ROA.
P mean cost
max-vol 10.4253
max-dist 10.1801
3-step nilpotent 10.3566
2-step nilpotent 10.3485
Table 5.3: Mean costs for Monte Carlo simulation using different policies
5.5.5 Heuristic Optimisation
Instead of solving directly for an inner approximation to the volume, a computationally simpler
problem can be solved that avoids the projection operation entirely. Given that the vector s
represents the amount of tightening applied to the nominal sets, minimising its two-norm may,
in some way, be expected to reduce the amount by which the ROA is shrunk by the constraint
tightening policy. Specifically, the optimisation problem is given by
Pheur(N) =

min
P
‖s‖2
subject to:
∃{Z(j) ≥ 0a×q}, {Z¯(j) ≥ 0a×r} :
s satisfies (5.17) and (5.18).
(5.70)
100
5.5. EXAMPLES
Applying this policy on different two-state systems shows that it gives the same region of
attraction as applying max-vol, suggesting that it may provide a simpler way of calculating
a volume-optimal policy. However, the heuristic produces a smaller volume for higher state
dimensions. Consider the three-state system
A =
1 1 0.50 1 1
0 0 1
 , B =
0.1670.5
1
 , C =

1 0 0
0 1 0
0 0 1
0 0 0
 , D =

0
0
0
1
 . (5.71)
Given a maximum prediction horizon length of N¯ = 4 steps, the control objective is to steer
the system to the target set
T1 =
{
x
∣∣∣∣ [−1 −1− 2]T ≤ x ≤ [1 1 2]T} , (5.72)
whilst satisfying the output constraints
Y =
{
y
∣∣∣∣ [−20 −5 −3 −3]T ≤ y ≤ [20 5 3 3]T} (5.73)
and being robust to the bounded disturbance contained within the set
W =
{
w
∣∣∣∣ [−0.1 −0.2 −0.2]T ≤ w ≤ [0.1 0.2 0.2]T} . (5.74)
Repeating the same analysis as in Subsection 5.5.1, Table 5.4 shows the comparison of ROA
volumes. In this case, the VH-ROA volumes are identical to the 4-step volumes, so only one
set of data is listed. Whilst the heuristic outperforms the other policies considered in this
scenario, it still results in a smaller ROA volume than max-vol. However, this means that it
can provide a reasonable approximation to max-vol whilst being computationally tractable
in high dimensions. Using a 1-norm or ∞-norm on s in (5.70) gives only marginally smaller
volumes than the 2-norm, so using these norms can quicken computation even more, given
that an LP rather than an SOCP needs to be solved. This technique is used in Chapter 6 for a
7-state system.
P vol
(⋃4
i=1R(i,P)
)
nominal 911.8229
max-vol 598.4078
heuristic 549.0186
max-dist 545.6656
3-step nilpotent 488.8198
Table 5.4: ROA volumes over different policies for 3-state system
101
5. OPTIMAL CONSTRAINT TIGHTENING POLICIES
5.6 Discussion
The examples in the preceding sections showed that max-vol performs well for two and three-
state systems. However, it is prudent to consider what may happen in higher dimensions. There
are both theoretical and computational issues which may affect performance and viability, and
these issues are addressed in this section. Potential solutions to these issues are also discussed.
5.6.1 Curse of Dimensionality
A potential theoretical problem arises with the inner-ellipsoidal approximation in higher
dimensions. As discussed in [Hay11], the ratio of the volume of a n-dimensional hypersphere
to its circumscribed n dimensional hypercube is given by
pi
n
2
2nΓ( n2 + 1)
, (5.75)
where Γ( · ) represents the gamma function. As n increases, this ratio becomes extremely small.
It can be reasonably expected that a similar phenomenon will be observed with polytopes of
certain shapes, so the inner ellipsoidal approximation will occupy only a small fraction of the
overall polytope volume, which tends to be concentrated near the vertices. One potential way
of addressing this issue is to use inner superellipsoids instead, with the same norm parameter
on each dimension. These are convex sets of the form,
E˜(c,Ω, p) = {x | x = c +Ωξ, ∀ ‖ξ‖p ≤ 1}, (5.76)
for centre c, shape matrix Ω  0 and norm parameter p. For computational tractability,
the restriction p ∈ {1, 2,∞} is imposed. Then, for p ∈ {1,∞}, the resulting shapes are
parallelotopes. In order to use this approximation, the definition of the matrix norm allows the
constraints (5.53) to be modified to
ΦˆTi c + ‖ΩΦˆi‖q ≤ ψˆi − Γis, (5.77)
where q = ∞ if p = 1 and q = 1 if p = ∞. Then, for certain polytope geometries, values of p
other than 2 may result in a better approximation to the true volume in higher dimensions.
Experimentation with the examples covered in the previous section for p ∈ {1,∞} gives the
same ROA as using the ellipsoid.
5.6.2 Computational Complexity
For higher dimensional problems, both the projection operation described in Theorem 5.5
and the optimisation problem (5.54) can pose a very high computational burden. For the
projection operation, an algorithm such as ESP whose complexity grows linearly with facet
count becomes necessary. In terms of reducing the complexity of the optimisation problem,
102
5.6. DISCUSSION
the heuristic optimisation approach proposed in Subsection 5.5.5 results in a computationally
simpler problem whilst still improving the region of attraction.
Another strategy to simplify the optimisation problem is to restrict the degrees of freedom in
which the ellipsoid can vary, by a priori specifying the directions in which the semiaxes will
point. This is possible due to the fact that Ω is positive definite, so it can be diagonalised to
give
Ω = TΛTT, (5.78)
where T is an orthogonal matrix and Λ is a diagonal matrix with the eigenvalues, which
correspond to the semiaxis lengths, on the diagonal. These eigenvalues will all be real and
positive. By specifying the matrix T, the search space is reduced. By also specifying a ratio
between the semiaxis lengths, the problem will collapse into a linear program. Using a
superellipsoid with p = 1 or p = ∞ with a fixed semiaxes ratio and directions is similar to the
approach taken in [BFT04], where inner-box approximations to polytopes are specified. Note
that specifying the orientation and aspect ratio of the ellipsoid a priori will require knowledge
of the typical ROA shape over varying tightening policies for a given system, which can be
difficult in higher dimensions.
5.6.3 Volume Calculation
An alternative to a convex problem that may be effective at higher dimensions is directly
applying nonlinear optimisation. By creating a function that returns the exact volume of a
polytope for a given tightening policy, nonlinear optimisation could be applied to find a locally-
optimal policy to maximise a fixed horizon ROA. To test this idea, projection and Delaunay
triangulation are used to define a volume function for the two dimensional problem defined
in Subsection 5.5.1. Applying nonlinear optimisation with different initial guesses does not
seem to find a better volume than that given by max-vol. However, not all starting guesses
result in the optimisation ultimately converging to the volume corresponding to max-vol
demonstrating the existence of local maxima. The issues with extending this approach to
higher dimensions are the projection operation and the exact volume computation for the
polytope, which can become computationally intensive, given that Delaunay triangulation
requires vertex enumeration as a preliminary step.
An alternative to exact volume measurement is estimation. Monte Carlo methods could be used,
which would allow direct estimation of the VH-ROA by counting the number of initial states
for which a feasible solution exists. This method still suffers from the curse of dimensionality
as well, as a very large number of random initial states will be needed for adequate coverage
in higher dimensions. Other methods include box-approximation algorithms [BFT04] and
entropy methods detailed in [BH10]. These noisy approximations could be used in conjunction
with gradient descent algorithms such as simultaneous-perturbation stochastic approximation
[Spa92] to locally optimise the policy for a fixed horizon ROA.
103
5. OPTIMAL CONSTRAINT TIGHTENING POLICIES
5.7 Conclusion
This chapter has developed a technique for choosing constraint-tightening policies that increase
the size of the region of attraction. By expressing the tightened sets as affine functions of the
constraint tightening policy and dual variables, a convex representation of the region of
attraction, parameterised by the tightening policy, is formed. A convex optimisation problem
is then posed to find the constraint tightening policy that maximises the volume of an inner-
ellipsoidal approximation to the ROA over fixed horizon lengths. Simulations with example
systems show this technique to enlarge even the full region of attraction over the variable
horizon, as well as permit vastly smaller prediction horizon lengths to be used when compared
to using nilpotent policies.
Future work will look at ways of formulating the optimisation problem for an inner approxima-
tion of the VH-ROA as a quasi-convex optimisation problem. It will further look at addressing
some of the issues that occur in higher state dimensions. Another avenue of investigation will
be the design of cost optimal policies. These extensions will potentially allow the optimisation
of a policy that provides a trade-off between ROA volume and mean cost over a given set of
starting states.
5.8 Acknowledgements
The work in this chapter has been submitted for review to the 50th IEEE Conference on
Decision and Control. The simulations in this chapter were carried out using MATLAB R2011b,
YALMIP [Löf04], the Multi-parametric Toolbox [KGB04], CPLEX [CPL09], SDPT3 [TTT99] and
the Invariant Set Toolbox [Ker00].
104
Chapter 6
Case Study - Surface Excavation
6.1 Introduction
Surface mining operations involve the recovery of valuable mineral ores from the earth through
removal of the covering material, known as overburden. In most cases, this requires the use of
high explosives in order to blast the overlying earth and ore to facilitate their removal. The
recovered ore and waste are separated and crushed for further processing.
The majority of surface mines utilise large-scale machinery for removing the blasted ore and
waste. These machines are presently controlled by skilled operators who require many years
of training in order to achieve proficiency. In particular, difficult digging conditions require
operators to have an intuitive understanding of material and bucket interactions. Control
strategies to be used in such conditions are only learnt through years of experience. This
experience results in the ability to achieve consistent bucket payloads, which is highly desirable
when considering the entire supply chain. These factors make the digging cycle of the machine,
where it engages material, an ideal candidate for automation.
MPC is a means through which such automation can be potentially achieved. Its ability to
handle constraints and optimise an objective function for performance suggest that it may
be useful for automating a surface mining machine. In this application, there is no notion of
stability once the bucket engages the material: the overall objective is to fill the bucket to a
target payload as quickly as possible and then exit the material. VH-MPC is therefore lends
itself to this control problem.
In order to demonstrate the applicability of VH-MPC to digging, the rope shovel is chosen as
a candidate machine for analysis. It is widely used in surface mining operations. As shown
in Figure 6.1, rope shovels have two degrees of freedom: the translational movement of the
crowd arm and the length of the hoist rope. The crowd arm pivots about the crowd pin. A
combination of hoist and crowd movements generate digging trajectories through material, as
shown in Figure 6.2
This chapter develops a simplified dynamic model of a representative rope-shovel with a
60-100 tonne nominal bucket capacity. A material model based on Mohr-Coulomb theory is
105
6. CASE STUDY - SURFACE EXCAVATION
θ
r
Crowd Pin
Boom Sheave
Bucket Attachment
Crowd Arm
Hoist Rope
Hoist Attachment
Figure 6.1: The rope shovel
Figure 6.2: Rope shovel excavation trajectory [P&03]
106
6.2. MECHANISM MODEL
also developed, simulating excavation in oil sand. The combined model of the shovel and
material goes further than [FHAO05] by including the inertial effects of a varying bucket
mass and allowing the material force to influence the bucket’s trajectory. Based on the overall
nonlinear model of the machine and material, an inner pre-stabilising controller is designed,
following which, a linear discrete-time model is identified. This model is used for designing
a VH-MPC controller for trajectory generation. Subsequently, simulations are performed,
which incorporate the theory on feasible contingencies, move blocking and optimal constraint
tightening developed in the previous chapters. These simulations show the generation of
realistic trajectories, motivating new applications of VH-MPC for autonomy that go beyond
the realm of vehicle manœuvring.
6.1.1 Nomenclature
Vectors denoting actual physical quantities such as positions, velocities, forces and moments
will be typeset in boldface. When the symbol appears without boldface, it denotes the mag-
nitude of the vector in question. Boldface will not be used for vectors in generalised coordinates
(where the components may have different physical units). Standard robotics notation will be
used to denote vectors in a given coordinate frame, so the leading subscript in nzm indicates the
vector zm in coordinate frame n. Position vectors will also be expressed in polar and Cartesian
form using the notation
nzm = rm
cos (
nθm)
sin (nθm)
0
 =

nxm
nym
0
 . (6.1)
Note that a different font is used to indicate cartesian components to distinguish them from
state and output variables. Although only two-dimensional motion is considered in this model,
the third vector component is required for the calculation of moments using cross-products.
The standard Cartesian unit vectors i, j and k will be used as needed.
6.2 Mechanism Model
A simplified rigid-body model of the shovel mechanism combined with bucket-material
interaction is developed for simulation purposes, based on reasonably realistic parameter
values from [FL07]. The kinematic and dynamic parameters are representative of a typical
medium capacity rope shovel. Two coordinates are selected as shown in Figure 6.1: crowd
arm translation (r) and crowd arm rotation (θ). These coordinates are kinematically related
to the hoist rope length. The various lengths and angles are shown in Figure 6.3, and the
corresponding parameter values are listed in Table 6.1. Also pictured are two systems of
coordinates used for specifying vectors: the base frame (0x,0y) fixed to the machine body and
(1x,1y) rotating with the crowd arm. The base frame is positioned 6.6 m above ground level.
Figure 6.4 illustrates the simplified mechanism model as well as the coordinates r and θ. The
107
6. CASE STUDY - SURFACE EXCAVATION
crowd arm and bucket are treated as a rigid body (combined mass M, inertia about combined
centre of mass J), with a sliding contact at the crowd pin approximating the rack-and-pinion on
an actual shovel. In addition, the weight of the material in the bucket, mp, is treated as a point
mass located at the bucket’s centre of mass. Estimated dynamic parameters, masses, centres of
mass, inertias and coefficients of friction are listed in Table 6.2.
0z2
0z4
0z1
0z3
1z6
1z7
1z8
1z50y
0x
1y
1x
Figure 6.3: Kinematic vectors
6.2.1 Kinematic Equations
For convenience, vectors in the base coordinate frame will have their leading superscripts
omitted. In terms of the crowd length r, pictured in 6.4, the vector 1z6 is defined as
1z6 =
[
r 0 0
]T
. (6.2)
Using the rotation matrix
Rθ =
cos θ − sin θ 0sin θ cos θ 0
0 0 1
 . (6.3)
which is a function of the crowd angle, this vector can be transformed into the base frame to
give
z6 = Rθ1z6 =
[
r cos θ r sin θ 0
]T
, (6.4)
The two degrees of freedom will be used as generalised coordinates for deriving equations
of motion in Subsection 6.2.2. Before doing this, however, kinematic equations are required
108
6.2. MECHANISM MODEL
M, J
mp
r
θ
Figure 6.4: Mechanism model showing generalised coordinates
Parameter Value Description
Lc 13.19 m Crowd arm length
r1 20.55 m Machine origin to boom sheave centre
r2 1.53 m Boom sheave radius
r4 7.39 m Machine origin to crowd pin
r5 1.47 m Crowd pivot radius about pin
r7 2.87 m Bucket attachment to hoist attachment
r8 7.64 m Bucket attachment to bucket teeth
rccm 6.7570 m Crowd pin to crowd’s centre of mass when r = Lc
rkcm 3.2000 m Bucket pin to bucket’s centre of mass
0θ1 41.74◦ Angle of vector z1 in base frame
0θ2 −22.50◦ Angle of vector z2 in base frame
0θ4 39.92◦ Angle of vector z4 in base frame
1θ5 90.00◦ Angle of vector z5 in crowd pin frame
1θ6 0.00◦ Angle of crowd arm in crowd pin frame
1θ7 39.42◦ Angle of vector z7 in crowd pin frame
1θ8 16.20◦ Angle of vector z8 in crowd pin frame
1θccm 12.57◦
Angle of vector from crowd pin to crowd centre of mass in crowd
pin frame
1θkcm −26.57◦ Angle of vector from bucket attachment to bucket centre or massin crowd pin frame
Table 6.1: Kinematic parameters
109
6. CASE STUDY - SURFACE EXCAVATION
Parameter Value Description
mc 65800 kg Mass of crowd arm
mk 24500 kg Mass of bucket
Jc 774900 kg m2 Crowd arm inertia about its mass centre
Jk 64580 kg m2 Bucket inertia about its mass centre
cr 1× 105 N s m−1 Viscous friction coefficient for crowd translation
cθ 1× 106 N m s Viscous friction coefficient for crowd rotation
g 9.81 m s−2 Acceleration due to gravity
Table 6.2: Dynamic parameters
to relate the hoist length and velocity to the generalised coordinates. This is due to the fact
that the crowd angle is in reality a function of both the crowd length and hoist rope length. In
addition, equations relating the generalised coordinates to the bucket tip position and velocity
are needed for trajectory information.
Closing the vector loop in Figure 6.3 gives the hoist vector relationship
z3 = z1 + z2 − z4 − Rθ
(
1z5 + 1z6 + 1z7
)
. (6.5)
The right hand side of (6.5) is a function of r and θ. The hoist length is then given by z3. From
the expression for z3, it is possible to relate the hoist rope velocity to r˙ and θ˙. By defining the
vector of generalised coordinates
q =
[
r θ
]T
, (6.6)
the time derivative of the hoist rope vector is calculated as
z˙3 = ∇qz3q˙
=
− cos θ r sin θ + r5 sin
(
θ + 1θ5
)
+ r7 sin
(
θ + 1θ7
)
− sin θ −r cos θ − r5 cos
(
θ + 1θ5
)− r7 cos (θ + 1θ7)
0 0
 q˙. (6.7)
In this expression, ∇qz3 is the vector gradient, or Jacobian, of the vector z3 with respect to the
generalised coordinates. Given that z3 =
√
z3 · z3, the chain rule can be applied to calculate the
rate of hoist extension or retraction
dz3
dt
=
z3
z3
· z˙3
= zˆ3 · z˙3. (6.8)
This is simply the component of z˙3 in the direction of the hoist unit vector zˆ3, which is expected
intuitively.
The position vector of the combined centre of mass 1zcm of the crowd arm and bucket, in the
110
6.2. MECHANISM MODEL
crowd pin frame is given by
1zcm =
1
M
(
mc1zccm +
[
r− Lc 0 0
]T
+ mk1zkcm + 1z5 + 1z6
)
, (6.9)
Note that this centre of mass moves depending on the crowd extension or retraction length.
The r− Lc term accounts for the fact that the position of the crowd centre of mass is defined
with the crowd fully extended. In terms of the parameters quantified in Table 6.1, the position
of the centre of mass in the base frame is given by the expression
z¯ = z4 + Rθ1zcm. (6.10)
Its velocity is then given by
˙¯z = ∇qz¯q˙, (6.11)
where ∇qz¯ is the Jacobian of the centre of mass with respect to the generalised coordinate
vector. Additionally, the velocity of the payload, assumed to be located at the bucket’s centre
of mass, is required. The position of the payload centre of mass in the base frame is then
z¯p = z4 + Rθ
(
1z5 + 1z6 + 1zkcm
)
. (6.12)
As with the equivalent centre of mass, the velocity of the payload is given by
˙¯zp = ∇qz¯pq˙. (6.13)
6.2.2 Equations of Motion
As a first step towards deriving the equations of motion for the mechanism, it is necessary to
find the combined inertia for the rigidly connected bucket and crowd-arm for in-plane rotation.
As this is independent of the coordinate frame, assume that the crowd arm is fully extended.
With this assumption, denote the resulting centre of mass position in the crowd pin frame as
1z′cm. This vector is constant in this frame. Then, the parallel axis theorem gives the equivalent
inertia of the combined crowd pin and bucket system about this centre of mass as
J = Jc + mc
((
1x′cm − 1xccm
)2
+
(
1y′cm − 1yccm
)2)
+ Jk + mk
((
1x′cm − Lc − 1xkcm
)2
+
(
1y′cm − r5 − 1ykcm
)2)
. (6.14)
Using the kinematic and dynamic parameters, the total kinetic energy of the system, denoted T
is calculated as
T =
1
2
M ˙¯z · ˙¯z+ 1
2
mp ˙¯zp · ˙¯zp + 12 Jθ˙
2. (6.15)
The potential energy, denoted V, depends only on the y components of the mechanism and
payload centres of mass in the base frame, so
V = Mgj · z¯+ mpgj · z¯p. (6.16)
111
6. CASE STUDY - SURFACE EXCAVATION
From the definition of the Lagrangian
L = T −V, (6.17)
the Euler-Lagrange differential equation is applied to each generalised coordinate to give
d
dt
(
∂L
∂r˙
)
− ∂L
∂r
= Fr − Fr, f (6.18)
d
dt
(
∂L
∂θ˙
)
− ∂L
∂θ
= τθ − τθ, f . (6.19)
The generalised force and moment Fr and τθ , as well as the dissipative force terms Fr, f and τθ, f
are quantified in terms of applied forces in Subsection 6.2.3, allowing the equations of motion
to be be written in the form [SHV05]
D(q, mp)q¨ + C(q, q˙, mp, m˙p) + G(q, mp) +Λq˙ = Φ, (6.20)
where D( · , · ) is the inertia matrix, C( · , · , · , · ) is a vector of Coriolis terms and G( · , · ) is a
vector of gravity terms. The dissipative forces are contained in the vector
Λq˙ =
[
Fr, f τθ, f
]T
(6.21)
and the applied external forces are denoted by the vector
Φ =
[
Fr τθ
]T
. (6.22)
These applied forces are transformations of the hoist rope and crowd arm actuator forces, as
well as the material resistance force Ftip, which is assumed to act at the bucket tip. As the
dynamic equations (6.20) depend on the varying payload mp, expressions for the matrices in
terms of the kinematic and dynamic parameters are rather long. Hence, they are presented in
Appendix B, having been derived using MATLAB’s symbolic toolbox. The expressions for the
external forces and frictional forces are derived in the following subsection.
6.2.3 External & Dissipative Forces
The equations of motion derived above require external forces to be expressed as generalised
forces acting on the coordinates r and θ. As frictional forces have not yet been quantified, they
will also be introduced here. Define the base-frame unit vector in the crowd extension direction
rˆ =
[
cos θ sin θ 0
]T
(6.23)
Then, the crowd extension force Fc, hoist pull force Fh, and tip force Ftip give the generalised
crowd force
Fr = Fc + rˆ ·
(
Fhzˆ3 + Ftip
)
. (6.24)
112
6.3. MATERIAL INTERACTION MODEL
Varignon’s theorem [BF08] is used to calculate the generalised moment about the crowd pin,
which is given by
τθ =
[
(z5 + z6 + z8)× Ftip + (z5 + z6 + z7)× Fhzˆ3 + τtip
] · k, (6.25)
where τtip is any moment applied to the bucket teeth. In the material model presented in the
next section, this is assumed to be zero, as the line of action of the force is assumed to pass
through the bucket teeth. The viscous friction term acting against crowd motion is given by
Fr, f = cr r˙. (6.26)
To quantify the dissipative moment it is necessary to sum both the viscous moment on the
crowd pin itself and the moment induced by the viscous friction acting on the crowd arm. The
resulting moment is given by
τθ, f = cθ θ˙ + crr5r˙. (6.27)
Hence, the dissipative force matrix is given by
Λ =
[
cr 0
r5cr cθ
]
(6.28)
6.3 Material Interaction Model
In general, analytically modelling bucket interaction behaviour is not possible. However,
homogeneous materials like oil sand can be approximated through Mohr-Coulomb theory,
based on static material failure. The Fundamental Equation of Earthmoving Mechanics (FEE)
[Ree64; LSC98; CS99] uses this theory to estimate the force required for a planar tool to
breach the shear failure threshold of the material, as a function of tool and terrain orientation.
Figure 6.5 illustrates the required parameters. Numerical values for the material parameters
can be found in Table 6.3 [Cla57; Col07]. The remainder are as follows: Lt is the length of the
tool extending under the terrain surface, d is the depth of the tool tip below the surface, L f is
the length of the failure plane below the surface, Q is the surcharge, or displaced soil that rests
on the wedge, φ is the material angle of repose, R is the resistance force of the material exerted
on the wedge. The material in the shaded wedge amounts to the payload i.e. the current bucket
fill. Its volume is denoted Vs.
In terms of the material and geometrical parameters, the force required to induce shear failure
is given by
Rs =
(
d2wγgNw + cwdNc +VsγgNq
)
hˆ, (6.29)
where w is the width of the tool and g is the acceleration due to gravity. Nw, Nc and Nq are
dimensionless parameters defined by the equations
Nw =
(cot β− tan α)(cos α+ sin α cot(β+ φ))
2[cos(ρ+ δ) + sin(ρ+ δ) cot(β+ φ)]
(6.30)
113
6. CASE STUDY - SURFACE EXCAVATION
Q
ρ
δ
α
β
ϕd
cLf
caLt
R
Rs
bucket blade
terrain surface
failure plane
Figure 6.5: Static balance of material failure forces [CS99]
Nc =
1+ cot β cot(β+ φ)
cos(ρ+ δ) + sin(ρ+ δ) cot(β+ φ)
(6.31)
Nq =
cos α+ sin α cot(β+ φ)
cos(ρ+ δ) + sin(ρ+ δ) cot(β+ φ)
. (6.32)
The direction of Rs is given by the unit vector
hˆ =
[
sin(ρ− α+ δ) cos(ρ− α+ δ) 0
]T
, (6.33)
and its line of action is, for the purposes of this model, assumed to pass through the bucket
teeth. Certain geometric configurations of the bucket can result in Rs becoming negative,
implying a material force pulling the bucket deeper horizontally, which is not observed in
reality. Hence, under such conditions, Rs is taken to be zero. Also note that for simplicity, a
linear terrain face profile defined by the equation[
1.66 −1
]
ztip = 27.018 m (6.34)
is used for simulation, where ztip ∈ R2 defines the (x, y) bucket tip position in frame 0. The fill
volume Vs is estimated from the bucket width and area subtended by the trajectory and face,
assuming the bucket were to pull out of the face vertically from its current pose as in Figure 6.5.
Given that the FEE only describes static failure conditions, it can be considered as a threshold
on the force that the material is able to resist until failure. Assuming that the sum of the forces
114
6.4. SYSTEM IDENTIFICATION
Parameter Value Description
α 59◦ Terrain surface (face) angle
β 35◦ Failure surface angle (assumed)
δ 37◦ Material-tool friction angle
φ 50◦ Angle of repose
γ 2000 kg m−3 Bulk density
c 0 Pa Cohesion
ca 0 Pa Material-tool adhesion
Table 6.3: Material properties
being exerted by the actuators and gravity at the bucket teeth is given by the vector Ft, the
resulting force being exerted by the soil at the bucket teeth is given by
Ftip = −min
(
Ft · hˆ, Rs
)
hˆ, (6.35)
where the force is exerted opposite to the failure direction. This transforms through (6.24) and
(6.25) to give a complete model of the machine digging through the material.
6.3.1 Implementation
In order to implement the dynamic model described by (6.20), it is necessary to define five
state variables, namely r, θ, r˙, θ˙ and mp. Then, define the state vector χ as
χ(t) ≡
[
r(t) θ(t) r˙(t) θ˙(t) mp(t)
]T
. (6.36)
Equation (6.20) can therefore be rewritten as a nonlinear state space equation of the form
χ˙(t) = f (χ(t),Φ(t)), (6.37)
where f ( · , · ) is a nonlinear vector function. This is now in a suitable form for simulation with
standard ODE solvers.
6.4 System Identification
As a first step towards identifying a linear model for the shovel, a pre-stabilising linear
controller with integral action is designed for trajectory tracking, to help mitigate some of the
effects of the nonlinear dynamics. This linear controller is designed by sampling data from
the open-loop nonlinear system to identify a suitable linearisation, then applying LQR-like
methods. A linear model is subsequently identified for the closed loop nonlinear trajectory-
tracking system, which will be used for MPC.
115
6. CASE STUDY - SURFACE EXCAVATION
6.4.1 Pre-Stabilising Controller
Given that there are no open-loop equilibrium states that are of interest for digging control,
a forced equilibrium is calculated, for an initial machine pose at rest given by r0 = 6 m and
θ0 = −60◦, with the bucket tip positioned just within the material face. Perturbations are made
to the force and moment keeping the crowd arm in equilibrium, causing the machine to move
through the material. During the resulting motion, data for four states, namely r, θ, r˙ and θ˙,
are collected with a sampling frequency 100 Hz. This allows a discrete-time linearisation to be
found using least-squares techniques, with the state and input vectors given by
x(k) =
[
r˜(k) θ˜(k) r˙(k) θ˙(k)
]T
(6.38)
u(k) =
[
F˜c(k), F˜h(k)
]T
, (6.39)
where t = kTs. The tilde on the states and inputs indicates a trim value from equilibrium. As
the machine is at rest in the trim condition, the crowd velocity and angular velocity correspond
to the actual states.
Remark 6.1. The trim value for the crowd force has an important role to play in correcting an
inadequacy of the material model. Given that the model does not account for the supporting
force which would act on the base of the bucket, the crowd is assumed to be acted upon by the
equilibrium force at all times. This simplification results in more realistic crowd forces than
ignoring the supporting force entirely.
To achieve offset-free tracking, this linear model is augmented with two discrete-time integral
states. In terms of the hoist and crowd position estimates from equilibrium, the integrator has
dynamics
η(k + 1) = η(k) +
[
r˜(k) θ˜(k)
]T − [r˜d(k) θ˜d(k)]T , (6.40)
where η( · ) ∈ R2×1 is a vector of integral states, r˜d( · ) is the desired trimmed crowd length and
θ˜d( · ) is the desired trimmed crowd angle. A state-feedback controller is now designed using
an LQR-like approach, with some extra manual tuning of the integral gains. Assuming that the
desired (reference) trajectory inputs are zero, a stabilising gain matrix K is found such that the
feedback law
u(k) = −K
[
x(k)
η(k)
]
(6.41)
provides adequate tracking performance with sufficient damping on the crowd velocity and
angular velocity. The numerical value of the gain can be found in Appendix B. Figure 6.6
illustrates the controller’s performance when tracking a simple piecewise-linear reference.
Tuning to achieve a non-oscillatory response results in the r reference tracking delay. It is
inconsequential, however, as the next stage of system identification allows it to be predicted.
116
6.4. SYSTEM IDENTIFICATION
0 5 10 15 20 25 30 35 40
5.5
6
6.5
7
7.5
time (s)
r 
(m
)
0 5 10 15 20 25 30 35 40
−1.2
−1
−0.8
−0.6
−0.4
time (s)
θ
 
 
Actual
Setpoint
Figure 6.6: Performance of trajectory tracking controller
6.4.2 Identification of Trajectory Tracking System
Having stabilised the nonlinear system, a linear model is now identified for the closed-loop
system, taking the trajectory reference as its input. Using the MATLAB system identification
toolbox, a subspace method [VODM96] with prediction error minimisation is employed to
identify the model, where it is assumed that there is no input feedforward. This time, the
payload and integral states are also included, producing a linear model with seven states. A
better linearisation is obtained by taking the square root of the payload trajectory data before
identification, as this state varies quadratically with linear changes in r.
The outputs are chosen to feed forward all of the states, in addition to the Cartesian bucket
teeth position deviations from the equilibrium pose. Note that the payload and position of the
bucket teeth depend nonlinearly on r( · ) and θ( · ), but predictions are required for enforcing
constraints later on with MPC, so linear estimates are desirable. The purpose of including the
integral states in the linear model is to allow constraining of actuator forces, since they can be
inferred from −Kx.
The linearised closed-loop model adequately predicts most of the states for unseen inputs
when compared to the original system. In order to utilise it for trajectory generation, it is
downsampled to 1 Hz using a zero-order hold, giving a new discrete time system linear system
of the form
xˆ(k + 1) = Axˆ(k) + Buˆ(k) (6.42)
yˆ(k) = Cxˆ(k) + Duˆ(k), (6.43)
117
6. CASE STUDY - SURFACE EXCAVATION
where A ∈ R7×7, B ∈ R7×2, C ∈ R9×7 and D = 09×2. Index k is redefined accordingly to
account for the new sampling frequency. The input, state and output vectors are defined by
uˆ(k) =
[
r˜d(k) θ˜d(k)
]T
(6.44)
xˆ(k) =
[
r˜(k) θ˜(k) r˙(k) θ˙(k)
√
mp(k) ηT(k)
]T
(6.45)
yˆ(k) =
[
xˆT(k) z˜Ttip(k)
]T
. (6.46)
Note that the output vector contains all of the states along with an estimate of the bucket tip
position from the trim condition, namely z˜tip(k).
The performance of the downsampled model is illustrated in Figure 6.7, for the crowd length
and angle. The performance in tracking all other signals and outputs is similarly close, apart
from the payload state, which is pictured in Figure 6.8.
0 5 10 15 20 25 30 35 40
5.5
6
6.5
7
7.5
time (s)
r 
(m
)
0 5 10 15 20 25 30 35 40
−1.2
−1
−0.8
−0.6
−0.4
time (s)
θ
 
 
Actual
Linearisation
Figure 6.7: Performance of identified linear model for crowd length and angle
6.5 VH-MPC Controller Design
The result of the preceding identification is a linear model that estimates the true state and
bucket tip trajectory of the nonlinear system, given the desired hoist and crowd position
signals as inputs. However, costing these inputs directly is not useful. Instead, the model must
be modified to allow costing of the change in u( · ) from its previous value, by augmenting
118
6.5. VH-MPC CONTROLLER DESIGN
0 5 10 15 20 25 30 35 40
0
20
40
60
80
100
120
time (s)
Pa
ylo
ad
 (t)
 
 
Actual
Linearisation
Figure 6.8: Performance of identified linear model for payload
additional states [Mac02] to store u(k− 1). This gives the prediction model[
x(k + j + 1|k)
u(k + j)
]
=
[
A B
0 I
] [
x(k + j|k)
u(k + j− 1|k)
]
+
[
B
I
]
∆u(k + j|k) (6.47)
y(k + j|k) =
[
C D
] [ x(k + j|k)
u(k + j− 1|k)
]
+ D∆u(k + j|k), (6.48)
where ∆u(k) = u(k)− u(k− 1). A maximum prediction horizon length N¯ = 15 steps is chosen,
as a typical digging cycle time for this material and machine is approximately 10-15 seconds.
6.5.1 Constraints
Constraints are now defined on the various system parameters. The critical constraints are
given by
|r˜| ≤ 5 m (6.49)
|r˜d| ≤ 4 m (6.50)
|∆r˜d| ≤ 0.2 m (6.51)
0 ≤ ∆θd ≤ 0.05 (6.52)
0 ≤ mp ≤ 100 t (6.53)∣∣F˜c∣∣ ≤ 5× 105 N (6.54)
119
6. CASE STUDY - SURFACE EXCAVATION
∣∣F˜h∣∣ ≤ 1.4× 106 N (6.55)[
−1 m 6 m
]T ≤ z˜tip ≤ [−2 m 15 m]T . (6.56)
The constraint on the crowd position reference rd is tighter than the actual constraint on
the crowd position r, to add an extra level of protection against violating this constraint.
Constraints on the change in crowd position and angle references ensure smooth trajectories.
The change in angle reference is constrained to be greater than zero, as during a standard
digging cycle, the hoist should only retract. The payload is limited to a maximum value of
100 t to avoid overloading the machine, and can never be negative. The crowd and hoist force
constraints ensure that these forces never reach the stall thresholds for those motors. Finally,
the tip position of the bucket is constrained to be within a reasonable region of the state space.
This final constraint is the reason that no bound constraints are needed on the angle of the
crowd arm. Note that, as with the states, the tilde on the tip position vector indicates a trimmed
value.
For dig completion, there are two terminal constraints imposed, in addition to all of the state
constraints in (6.49) - (6.56) being satisfied, namely
−30◦ ≤ θ ≤ −25◦ (6.57)[
1.66 −1
]
z˜tip ≤ −0.32 m. (6.58)
The former constraint implies that the crowd arm is at a suitably high angle for dig completion
and the latter ensures that the bucket tip has emerged from the face. Given that the linear
payload estimate is not highly accurate, no terminal payload constraint is imposed. However,
soft constraints are imposed through the cost function, as explained below.
As with the inverted pendulum example in Chapter 5, some constraint tightening is applied to
help mitigate the effect of linearisation error. Assuming appropriate units on each channel, a
disturbance vector of the form
W =
{
w
∣∣∣∣ |w| ≤ 1× 10−4 [1 1 1 1 1× 104 1 1 1 1]T} , (6.59)
is introduced, which applies a large amount of tightening to the payload and a smaller margin
of tightening to the remaining states. The payload disturbance is made as large as possible
whilst still ensuring non-empty tightened sets. This includes the tightened sets for feasible
contingencies and move blocking, which will be presented in Subsections 6.6.2 and 6.6.3. Given
the high state dimension and long prediction horizon, the heuristic optimisation approach
presented in Section 5.5, using a 1-norm cost, is used to calculate the policy. The region of
attraction is not critical for this scenario, as a fixed starting state has been assumed.
120
6.6. SIMULATION
6.5.2 Cost Function
The cost function is designed to weight the completion time with commanded input changes
and state-based operating costs, which are effectively soft constraints. The effect is to “shape”
the digging trajectory to encourage sufficient penetration depth of the material and adequate
payload acquisition. The cost function has the form
N(k)−1
∑
j=0
1+ 0.01 ‖u(k + j|k)‖1 + 0.8 ‖r˙(k + j|k)‖1 + 0.3
∥∥θ˙(k + j|k)∥∥1
+ 0.1
∥∥∥[1.66 −1] z˜tip(k + j|k)− 1.29∥∥∥
1
+ 0.01
∥∥∥√mp(k + j|k)−√65× 103∥∥∥
1
, (6.60)
which can be written in the form (3.40) by appropriately choosing the matrix Θ and output
reference yr. Notice the square root of the payload appearing in the last term of the cost
function, which is a result of having identified a system based on this quantity rather than the
payload itself. The cost function has the following effects:
• penalise time to completion
• penalise the value of the change in crowd position and angle reference
• penalise high crowd translational and rotational velocities to encourage smooth trajector-
ies
• encourage the face to be penetrated to a depth of 0.67 m.
• encourage a payload close to 65 t.
6.6 Simulation
Using the controller designed in the previous section, the system is simulated, starting with all
states of the prediction model initialised to zero (i.e. at their trim values for the actual model).
The first set of simulations consider a normal digging scenario. The next two sets consider the
implementation of a feasible contingency and move blocking respectively.
6.6.1 Normal Operation
Figure 6.9 shows the predicted tip trajectory from the linear model plotted with the actual
trajectory obtained when using the VH-MPC controller on the stabilised nonlinear model. It
can be seen that a realistic trajectory is generated. The trajectories for the crowd position and
angle are shown in Figure 6.10, together with the reference signal generated by the VH-MPC
controller. The payload acquisition, shown in Figure 6.11 is slightly less than the target value,
121
6. CASE STUDY - SURFACE EXCAVATION
but this is partly due to the error in estimating this state in the linear model. As shown in
Figure 6.12, the actuator forces respect the constraints imposed on them. The components of
the cutting force are displayed in Figure 6.13. For the actuator forces and the cutting force, the
values are shown only for the time whilst the bucket is in the material.
14 16 18 20 22
−7
−6
−5
−4
−3
−2
−1
0
1
2
3
x (m)
y 
(m
)
 
 
Trajectory
Prediction
Face
Figure 6.9: Actual trajectory and prediction from linear model
6.6.2 Stall Contingency
A contingency scenario is now implemented on the shovel, to handle a crowd stalling condition.
This could be caused by a large particle within the material, or a number of small particles
becoming lodged between the bucket teeth. Under such conditions, the crowd is unable to
extend and can only retract. Such stalls are typically handled by retracting the crowd whilst
still retracting the hoist. In the case of a severe stall, the bucket is to be withdrawn from the
face. The contingency therefore constrains the crowd reference to be non-negative and ensures
that a pull-back trajectory of at most 3 s exists to the face of the material.
Using the theory developed in Chapter 3, the same disturbance set (6.59) and optimised CT
policy is used for the primary and contingency constraints, which now allow a slightly higher
maximum magnitude for the change in crowd and hoist reference (0.3 m and 0.7 respectively).
Figure 6.14 illustrates the effect on the resulting trajectory under normal operation. As expected,
the shovel now digs closer to the face to ensure the existence of the pullout trajectory. The
actual pullout trajectory is illustrated in green, assuming a crowd stall at 6 s. Figure 6.15 shows
that the bucket fill is slightly lower when the contingency is required to be available.
122
6.6. SIMULATION
0 2 4 6 8 10 12
5
5.5
6
6.5
7
time (s)
r 
(m
)
0 2 4 6 8 10 12
−1.2
−1
−0.8
−0.6
−0.4
time (s)
θ
 
 
Trajectory
Reference
Figure 6.10: Crowd position and angle compared to reference
0 2 4 6 8 10 12
0
10
20
30
40
50
60
70
time (s)
Pa
ylo
ad
 (t)
Figure 6.11: Payload acquisition
123
6. CASE STUDY - SURFACE EXCAVATION
0 2 4 6 8 10
0
1
2
3
4
5
x 105
time (s)
Cr
ow
d 
Fo
rc
e 
(N
)
0 2 4 6 8 10
0
5
10
15
x 105
time (s)
H
oi
st
 F
or
ce
 (N
)
Figure 6.12: Hoist and crowd forces from the trim point
0 2 4 6 8 10
−4.5
−4
−3.5
−3
−2.5
−2
−1.5
−1
−0.5
0
x 105
time (s)
R
es
is
ta
nc
e 
Fo
rc
e 
(N
)
 
 
x component
y component
Figure 6.13: Cutting resistance force
124
6.6. SIMULATION
14 16 18 20 22
−7
−6
−5
−4
−3
−2
−1
0
1
2
3
x (m)
y 
(m
)
 
 
Normal Operation
With Stall Contingency
Contingency Activated
Face
Figure 6.14: Trajectory comparison with crowd stall contingency
0 2 4 6 8 10
0
10
20
30
40
50
60
70
time (s)
Pa
ylo
ad
 (t)
 
 
Normal Operation
With Stall Contingency
Figure 6.15: Payload comparison with crowd stall contingency
125
6. CASE STUDY - SURFACE EXCAVATION
6.6.3 Move Blocking
It is interesting to investigate the effect of move blocking on the closed loop trajectories. Using
the input and state condensed partitioned horizon strategy formulated in Chapter 4, with initial
sub-horizon length 5, terminal sub-horizon length 1 and medial sub-horizon block length 9, the
normal operation scenario is considered for simulation. The resulting trajectories are shown
in Figure 6.16 and the payloads are compared in Figure 6.17. In this case, a reduction in the
number of input variables by over 50% produces only a marginal reduction in the resulting
payload, but a longer completion time. The maximum solution time per timestep with move
blocking is 0.1814 s (serial computation), whereas without blocking it is 0.4437 s; this is a
reduction of almost 60% with blocking.
6.7 Conclusion
This chapter has demonstrated how VH-MPC can be applied to trajectory generation for surface
excavation using a rope shovel. It has also shown the application of the theory developed in
the previous chapters, namely feasible contingencies, move blocking and optimal constraint
tightening, to a higher-dimensional model with a longer prediction horizon. Reasonably
realistic trajectories have been obtained, confirming that VH-MPC has applicability beyond
vehicle path-planning problems.
6.8 Acknowledgements
The work in this chapter is based upon a paper presented by the author at the 49th IEEE
Conference on Decision and Control [SM10]. The simulations in this chapter were carried
out using MATLAB R2011b, YALMIP [Löf04], the Multi-parametric Toolbox [KGB04], CPLEX
[CPL09] and the Invariant Set Toolbox [Ker00].
126
6.8. ACKNOWLEDGEMENTS
14 16 18 20 22
−7
−6
−5
−4
−3
−2
−1
0
1
2
3
x (m)
y 
(m
)
 
 
Normal
Move Blocked
Face
Figure 6.16: Trajectory comparison with move blocking
0 5 10 15
0
10
20
30
40
50
60
70
time (s)
Pa
ylo
ad
 (t)
 
 
Normal
Move Blocked
Figure 6.17: Payload comparison with move blocking
127
128
Chapter 7
Conclusion
This thesis has developed new techniques and algorithms for robustness and optimality
in variable horizon model predictive control. These include: the development of control
algorithms robust to changes in system dynamics, constraints and objectives; an algorithm to
reduce computational complexity whilst maintaining robustness guarantees; and a technique
to optimise a controller’s region of attraction when robustness guarantees are required. This
chapter summarises these key contributions and motivates directions for future work.
7.1 Contributions
This thesis makes four main contributions to the body of academic knowledge. These contribu-
tions are summarised in the subsections below.
7.1.1 Feasible Contingencies
To address control scenarios where system dynamics, constraints and control objectives can
change, the notion of a controller with feasible contingencies has been introduced for comple-
tion problems. The key feature of this work is that, although the effect of the contingency is
known, the activation time is not known a priori. Criteria for nominal and robust contingency
availability at all times have been proposed. Sufficient conditions have then been used to design
a VH-MPC controller to fulfil these criteria, whilst maintaining robust recursive feasibility
and finite time completion guarantees. A double constraint tightening approach has been
formulated to ensure robustness before and after contingency activation.
Various contingency scenarios have then been presented, to encompass a wide gamut of po-
tential applications. These scenarios include multiple, costed, prioritised and state-dependent
contingencies. Two systems have then been used to showcase the application of the developed
theory. A simple vehicle example has been used to illustrate the feasible contingency idea on a
fault-tolerant control problem. In addition, a cart pendulum model has shown the applicability
of feasible contingencies to regulation problems.
129
7. CONCLUSION
7.1.2 Complexity-Reduced VH-MPC with Move Blocking
Given the computational complexity of VH-MPC optimisation problems over long time hori-
zons, new theory has been developed to apply input move blocking to completion problems,
to reduce computational complexity. A new notation has been proposed to describe the action
of move blocking; this has been subsequently applied to the design of a robust move-blocked
variable horizon controller. A modified technique for tightening constraints has then been
formulated, which preserves robust recursive feasibility and finite-time completion guarantees,
subject to certain conditions on the cost function.
As the new framework provides flexibility in the choice of blocking regime, a particular form
known as the partitioned horizon has been described, which enables fine-grained control action
at the beginning and end of a control manœuvre, whilst implementing move blocking in the
middle for complexity reduction. A simple example has then been used to illustrate the effect
of different block lengths when using this regime. From this example, it has been observed
that a reduction in decision variables, whilst reducing the computational burden as expected,
does not necessarily increase the cost.
7.1.3 Optimal Constraint Tightening Policies
A new algorithm for designing disturbance-feedback constraint tightening policies has been
presented, which improves the region of attraction for VH-MPC controllers. The algorithm
has been designed around the solution of a convex optimisation problem, affinely paramet-
erised in terms of the constraint tightening policy. By maximising the volume of an inner
ellipsoidal approximation to the true ROA, improvements in region of attraction size have been
demonstrated, in comparison to other popular approaches to designing CT policies. It has also
been shown how the algorithm proposed in this thesis can reduce the maximum prediction
horizon length. Methods have been proposed for dealing with high state-dimensions, where
the ellipsoid maximisation method may struggle computationally.
7.1.4 Case Study on Surface Excavation
A reasonably realistic model of a rope shovel used in surface excavation has been developed
to illustrate the efficacy of VH-MPC and the theoretical techniques developed in this thesis.
The model has applied VH-MPC to generate trajectories which are then tracked by an inner
pre-stabilising controller. The resulting realistic trajectories have motivated new applications
of VH-MPC for autonomy that go beyond the realm of vehicle path planning.
7.2 Future Research Directions
This section investigates future avenues of research stemming from the results presented in
this thesis.
130
7.2. FUTURE RESEARCH DIRECTIONS
7.2.1 Uncertain Contingencies
The work in this thesis assumes that the system dynamics after contingency activation are
well defined, apart from some bounded state disturbance. It is interesting to investigate to
what degree contingency availability can be guaranteed for uncertain systems, where the state
space matrices describing the dynamics after contingency activation are contained within some
polytope. This would be especially useful in fault-tolerant control scenarios, where the mode
of failure might not be known exactly.
7.2.2 Optimised & Multiplexed Move Blocking
The partitioned horizon has been proposed as an example blocking regime, but it has been
shown that the closed loop cost performance when applying the regime is not monotonically
related to the number of decision variables. This poses the problem of determining the closed-
loop cost-optimal blocking regime under disturbances. It may be difficult to determine this
in general, as the optimal regime is likely to be system and problem dependent. In this case,
some heuristics or rules of thumb for selecting blocking regimes for certain classes of problems
will be desirable.
Another avenue of research related to move blocking is the application of different regimes
on each input channel. More research is required to design a new constraint tightening
algorithm to maintain recursive feasibility and finite-time completion guarantees with this
scheme. Finding the optimal blocking regime for each channel also presents another research
opportunity.
7.2.3 Cost Optimal Constraint Tightening Policies
The algorithms presented for optimal constraint tightening in this thesis have focused on
improving the region of attraction. In some applications however, the expected closed loop
cost may be more important. A worthwhile avenue of investigation is policy optimisation for
some weighted sum of ROA size and expected closed loop cost. If the set of initial states is
small, then it would also be desirable to find the cost optimal policy where the ROA is just
large enough to contain this set.
7.2.4 Combinations of Strategies
Whilst the key contributions of feasible contingencies, move blocking and optimal constraint
tightening have been presented separately, further research will allow them to be combined.
Move blocking on the primary and contingency trajectories would reduce the complexity of
the variable horizon feasible contingency algorithm. Also, the optimal constraint tightening
algorithm could be modified to account for move blocking constraints, which would allow the
ROA volume to be maximised under blocking. The problem of finding the optimal combination
131
7. CONCLUSION
of blocking and constraint tightening is a future research direction that follows naturally from
this.
7.2.5 Shovel & Material Model Generalisation
In Chapter 6, a single linearisation has been obtained for the shovel and material system, which
assumed a particular face profile. To generalise this to arbitrary profiles, an algorithm is needed
to identify an appropriate linearisation online, profiling the current face prior to commencing
the digging cycle. The algorithm would also have to appropriately tune the cost function of the
VH-MPC controller, depending on the identified model.
To further generalise this algorithm to different materials, other modelling approaches will
be required. Whilst Mohr-Coulomb theory is suitable for homogeneous materials, it does not
adequately model blasted materials which can have particles on several different length scales.
To capture the behaviour of such materials, Discrete Element Modelling (DEM) can be used.
First proposed by [Cun71] to model the failure of brittle blocks of rock, DEM is an adaptation
of methods used widely in physics and chemistry to model liquid and gas behaviours [HF06].
Using a suitable software package, the dynamic model of the shovel could be interfaced to the
material model, and a linearisation identified for the combined system before applying MPC.
132
Bibliography
[BF08] A. Bedford and W. L. Fowler. Engineering Mechanics : Statics & Dynamics. 5th.
London: Pearson Education, 2008.
[BFT04] A. Bemporad, C. Filippi and F. D. Torrisi. “Inner and outer approximations of
polytopes using boxes”. In: Computational Geometry 27.2 (2004), pp. 151 –178.
[BH07] L. Breger and J. P. How. “Powered safe abort for autonomous rendezvous of
spacecraft”. In: AIAA Guidance, Navigation and Control Conference. 2007.
[BH08] L. Breger and J. How. “Safe trajectories for autonomous rendezvous of spacecraft”.
In: Journal of Guidance, Control and Dynamics 31.5 (2008), pp. 1478–1496.
[BH10] A. Barvinok and J. Hartigan. “Maximum entropy gaussian approximations for
the number of integer points and volumes of polytopes”. In: Advances in Applied
Mathematics 45.2 (2010), pp. 252 –289.
[BM99] A. Bemporad and M. Morari. “Control of systems integrating logic, dynamics,
and constraints”. In: Automatica 35.3 (1999), pp. 407 –427.
[BV04] S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press,
2004.
[Boy+94] S. Boyd, L. El Ghaoui, E. Feron and V. Balakrishnan. Linear Matrix Inequalities
in System and Control Theory. Philadelphia: Society for Industrial and Applied
Mathematics, 1994.
[CLM96] L. Chisci, A. Lombardi and E. Mosca. “Dual receding horizon control of con-
strained discrete-time systems”. In: European Journal of Control 2 (1996), pp. 278–
285.
[CPL09] CPLEX. IBM ILOG CPLEX V12.1 User’s Manual. 2009.
[CRZ01] L. Chisci, J. A. Rossiter and G. Zappa. “Systems with persistent disturbances:
predictive control with restricted constraints”. In: Automatica 37.7 (July 2001),
pp. 1019–1028.
[CS99] H. Cannon and S. Singh. “Models for automated earthmoving”. In: Experimental
Robotics VI, Lecture Notes in Control and Information Sciences. Springer Verlag, 1999,
pp. 183–192.
[Cag+07] R. Cagienard, P. Grieder, E. Kerrigan and M. Morari. “Move blocking strategies in
receding horizon control”. In: Journal of Process Control 17.6 (2007), pp. 563 –570.
133
BIBLIOGRAPHY
[Cam09] C Campbell. “Model Predictive Control for Stability Augmentation and Guidance
of a Flexible Micro Air Vehicle.” MSc Thesis. University of Cambridge, 2009.
[Car+08] J. M. Carson, B. Açikmes¸e, R. M. Murray and D. G. MacMynowski. “A robust
model predictive control algorithm with a reactive safety mode”. In: Proceedings
of the 17th IFAC World Congress. 2008.
[Cla57] K. Clark. Bulk Densities, Porosities and Liquid Saturations of Good Grade Athabasca
Oil Sands. Tech. rep. Research Council of Alberta, 1957.
[Col07] P. M. Collins. “Geomechanical effects on the SAGD process”. In: SPE Reservoir
Evaluation & Engineering 10.4 (2007), pp. 367–375.
[Com74] L Comtet. Advanced Combinatorics: The Art of Finite and Infinite Expansions, rev. enl.
ed. Dordrecht, Netherlands: Reidel, 1974, pp. 176–177.
[Cun71] P. Cundall. “Distinct element models of rock and soil structure”. In: Analytical and
Computational Methods in Engineering Rock Mechanics. Ed. by E. Brown. London:
Unwin Publishers, 1971, pp. 129–163.
[DE73] G. B. Dantzig and B. C. Eaves. “Fourier-motzkin elimination and its dual”. In:
Journal of Combinatorial Theory, Series A 14.3 (1973), pp. 288 –297.
[Del34] B. N. Delaunay. “Sur la sphère vide”. In: Bulletin of Academy of Sciences of the USSR
6 (1934), pp. 793–800.
[FHAO05] S. Frimpong, Y. Hu and K. Awuah-Offei. “Mechanics of cable shovel-formation
interactions in surface mining excavations”. In: Journal of Terramechanics 42.1
(2005), pp. 15 –33.
[FL07] S. Frimpong and Y. Li. “Stress loading of the cable shovel boom under in-situ
digging conditions”. In: Engineering Failure Analysis 14.4 (2007), pp. 702 –715.
[GI07] R. Gondhalekar and J. Imura. “Recursive feasibility guarantees in move-blocking
MPC”. In: Proceedings of the 46th IEEE Conference on Decision and Control. 2007,
pp. 1374–1379.
[GIK09] R. Gondhalekar, J. Imura and K. Kashima. “Controlled invariant feasibility –
a general approach to enforcing strong feasibility in MPC applied to move-
blocking”. In: Automatica 45.12 (2009), pp. 2869 –2875.
[GKM06] P. J. Goulart, E. C. Kerrigan and J. M. Maciejowski. “Optimization over state
feedback policies for robust control with constraints”. In: Automatica 42.4 (2006),
pp. 523 –533.
[GKR97] J. R. Gossner, B. Kouvaritakis and J. A. Rossiter. “Stable generalized predictive
control with constraints and bounded disturbances”. In: Automatica 33.4 (Apr.
1997), pp. 551–568.
[Glp] GNU Linear Programming Toolkit v4.47. 2010.
134
BIBLIOGRAPHY
[HF06] S. Hardy and E. Finch. “Discrete element modelling of the influence of cover
strength on basement-involved fault-propagation folding”. In: Tectonophysics
415.1-4 (2006), pp. 225–238.
[HM09] S. Heubach and T. Mansour. Combinatorics of Compositions and Words. CRC Press,
2009.
[Har10] E. N. Hartley. “Model Predictive Control for Spacecraft Rendezvous”. PhD thesis.
University of Cambridge, 2010.
[Hay11] B. Hayes. “An adventure in the nth dimension”. In: American Scientist 99 (2011),
pp. 442–446.
[JKM04] C. N. Jones, E. C. Kerrigan and J. M. Maciejowski. Equality Set Projection: A new
algorithm for the projection of polytopes in halfspace representation. Tech. rep. CUED/F-
INFENG/TR.463. Department of Engineering, University of Cambridge, 2004.
[KBM96] M. V. Kothare, V. Balakrishnan and M. Morari. “Robust constrained model predict-
ive control using linear matrix inequalities”. In: Automatica 32.10 (1996), pp. 1361
–1379.
[KG88] S. S. Keerthi and E. G. Gilbert. “Optimal infinite-horizon feedback laws for a
general class of constrained discrete-time systems: stability and moving-horizon
approximations”. In: Journal of Optimization Theory and Applications 57.2 (May
1988), pp. 265–293.
[KG95] I. Kolmanovsky and E. Gilbert. “Maximal output admissible sets for discrete-
time systems with disturbance inputs”. In: Proceedings of the American Control
Conference. Vol. 3. 1995, pp. 1995–1999.
[KGB04] M. Kvasnica, P. Grieder and M. Baotic´. Multi-Parametric Toolbox (MPT). 2004.
[KM04] E. C. Kerrigan and J. M. Maciejowski. “Feedback min-max model predictive
control using a single linear program: robust stability and the explicit solution”.
In: International Journal of Robust and Nonlinear Control 14.4 (2004), pp. 395–413.
[KRH07] Y. Kuwata, A. Richards and J. How. “Robust receding horizon control using
generalized constraint tightening”. In: Proceedings of the 26th American Control
Conference. 2007, pp. 4482–4487.
[Ker00] E. C. Kerrigan. “Robust Constraint Satisfaction: Invariant Sets and Predictive
Control”. PhD thesis. University of Cambridge, 2000.
[LSC98] O. Luengo, S. Singh and H. Cannon. “Modeling and identification of soil-tool
interaction in automated excavation”. In: Proceedings, IEEE/RSJ International Con-
ference on Intelligent Robots and Systems. Vol. 3. 1998, pp. 1900–1906.
[Lan+04] W. Langson, I. Chryssochoos, S. V. Rakovic and D. Q. Mayne. “Robust model
predictive control using tubes”. In: Automatica 40.1 (Jan. 2004), pp. 125–133.
[Löf04] J. Löfberg. “YALMIP : a toolbox for modeling and optimization in MATLAB”. In:
Proceedings of the CACSD Conference. Taipei, Taiwan, 2004, pp. 284 –289.
135
BIBLIOGRAPHY
[ML01] D. Mayne and W. Langson. “Robustifying model predictive control of constrained
linear systems”. In: Electronics Letters 37.23 (2001), pp. 1422 –1423.
[MM93] H. Michalska and D. Mayne. “Robust receding horizon control of constrained non-
linear systems”. In: IEEE Transactions on Automatic Control 38.11 (1993), pp. 1623–
1633.
[Mac02] J. M. Maciejowski. Predictive control with constraints. Essex, England: Prentice Hall,
2002.
[May+00] D. Mayne, J. Rawlings, C. Rao and P. Scokaert. “Constrained model predictive
control: stability and optimality”. In: Automatica 36.6 (2000), pp. 789 –814.
[May+11] D. Q. Mayne, E. C. Kerrigan, E. J. van Wyk and P. Falugi. “Tube-based robust non-
linear model predictive control”. In: International Journal of Robust and Nonlinear
Control 21.11 (2011), pp. 1341–1353.
[Old+09] F. Oldewurtel, R. Gondhalekar, C. Jones and M. Morari. “Blocking parameteriza-
tions for improving the computational tractability of affine disturbance feedback
MPC problems”. In: Proceedings of the 48th IEEE Conference on Decision and Control,
held jointly with the 28th Chinese Control Conference. 2009, pp. 7381–7386.
[RH03] A. Richards and J. How. “Model predictive control of vehicle maneuvers with
guaranteed completion time and robust feasibility”. In: Proceedings of the American
Control Conference. Vol. 5. 2003, pp. 4034–4040.
[RH06a] A. Richards and J. P. How. “Robust variable horizon model predictive control for
vehicle maneuvering”. In: International Journal of Robust & Nonlinear Control 16.7
(2006), pp. 333–351.
[RH06b] A. Richards and J. How. “Robust stable model predictive control with constraint
tightening”. In: Proceedings of the American Control Conference. 2006, pp. 1557–1562.
[RKC00] J. Rossiter, B. Kouvaritakis and M. Cannon. “Triple mode control in mpc”. In:
Proceedings of the American Control Conference. Vol. 6. 2000, pp. 3753–3757.
[RM09] J. B. Rawlings and D. Q. Mayne. Model Predictive Control: Theory and Design.
Madison, WI.: Nob Hill Publishing, 2009.
[Rak+06] S. Rakovic, E. Kerrigan, D. Mayne and J. Lygeros. “Reachability analysis of
discrete-time systems with disturbances”. In: IEEE Transactions on Automatic
Control 51.4 (2006), pp. 546 –561.
[Ree64] A. Reece. “The fundamental equation of earthmoving mechanics”. In: Proceedings
of the Institution of Mechanical Engineers. Vol. 179. 3F. 1964, pp. 16–22.
[Ric05a] A. Richards. “Robust model predictive control for time-varying systems”. In: 44th
IEEE Conference on Decision and Control and 2005 European Control Conference. 2005,
pp. 3747–3752.
136
BIBLIOGRAPHY
[Ric05b] A. G. Richards. “Robust Constrained Model Predictive Control”. PhD thesis.
Massachusetts Institute of Technology. Dept. of Aeronautics and Astronautics.,
2005.
[SHF04] T. Schouwenaars, J. How and E. Feron. “Receding horizon path planning with
implicit safety guarantees”. In: Proceedings of the 2004 American Control Conference.
Vol. 6. 2004, pp. 5576–5581.
[SHV05] M. W. Spong, S. Hutchinson and M. Vidyasagar. Robot Modelling and Control.
Wiley, 2005.
[SM10] R. C. Shekhar and J. M. Maciejowski. “Surface excavation with model predictive
control”. In: Proceedings of the 49th IEEE Conference on Decision and Control (CDC).
2010, pp. 5239–5244.
[SM11] R. C. Shekhar and J. M. Maciejowski. “Robust predictive control with feasible
contingencies for fault tolerance”. In: Proceedings of the 18th IFAC World Congress.
2011, pp. 4666–4671.
[SM12] R. C. Shekhar and J. M. Maciejowski. “Robust variable horizon MPC with move
blocking”. In: Systems & Control Letters 61.4 (2012), pp. 587 –594.
[SM98] P. Scokaert and D. Mayne. “Min-max feedback model predictive control for
constrained linear systems”. In: IEEE Transactions on Automatic Control 43.8 (1998),
pp. 1136–1142.
[SMR99] P. Scokaert, D. Mayne and J. Rawlings. “Suboptimal model predictive control
(feasibility implies stability)”. In: IEEE Transactions on Automatic Control 44.3
(1999), pp. 648–654.
[Sch+01] T. Schouwenaars, E. Feron, B. de Moor and J. How. “Mixed integer programming
for multi-vehicle path planning”. In: Proceedings of the European Control Conference.
2001.
[Sch86] A. Schrijver. Theory of Linear and Integer Programming. Chichester, UK: Wiley-
Interscience, 1986.
[Spa92] J. Spall. “Multivariate stochastic approximation using a simultaneous perturb-
ation gradient approximation”. In: IEEE Transactions on Automatic Control 37.3
(1992), pp. 332 –341.
[TTT99] K. C. Toh, M. Todd and R. H. Tütüncü. “SDPT3 - a MATLAB software package
for semidefinite programming”. In: Optimization Methods and Software 11 (1999),
pp. 545–581.
[Tro09] P. A. Trodden. “Robust Distributed Control of Constrained Linear Systems”. PhD
thesis. University of Bristol, 2009.
[VODM96] P. Van Overschee and B. De Moor. Subspace Identification for Linear Systems: Theory,
Implementation, Applications. Kluwer Academic Publishers, 1996.
137
BIBLIOGRAPHY
[WB10] Y. Wang and S. Boyd. “Fast model predictive control using online optimization”.
In: IEEE Transactions on Control Systems Technology 18.2 (2010), pp. 267 –278.
[ZG03] Y. Zhang and L. Gao. “On numerical solution of the maximum volume ellipsoid
problem”. In: SIAM Journal on Optimization 14.1 (2003), pp. 53–76.
[P&03] P&H MinePro Services. Peak Performance Practices: Excavator Selection. 2003.
138
Appendix A
Constraint Tightening Addenda
A.1 Derivation of the Max-Dist Policy
The aim of this policy is to replace the disturbance set with one of the form
W = {w | ζw ≤ βθ}, (A.1)
for some variable scaling factor β > 0. Given the effect of the scaling factor β on s( · ) and σ( · )
in (5.17) and (5.18), the non-emptiness constraints (5.61) need to be reformulated as
∃y, x : (A.2)
Ey + ey
( q
C
i=1
‖Ei‖
)T
≤ f − βs(N − 1)
Gx + ex
( q
C
i=1
‖Gi‖
)T
≤ h− βσ(N).
(A.3)
Then, this implies that
∃y, x : (A.4)
1
β
Ey +
1
β
ey
( q
C
i=1
‖Ei‖
)T
≤ 1
β
f − s(N − 1)
1
β
Gx +
1
β
ex
( q
C
i=1
‖Gi‖
)T
≤ 1
β
h− σ(N).
(A.5)
To retain convexity, define the variables δ = 1/β, y˜ = y/β and x˜ = x/β to give the equivalent
constraints
∃y˜, x˜ : (A.6)
Ey˜ ≤ δ
(
f − ey
( q
C
i=1
‖Ei‖
)T)
− s(N − 1)
Gx˜ ≤ δ
(
h− ex
( q
C
i=1
‖Gi‖
)T)
− σ(N).
(A.7)
139
A. CONSTRAINT TIGHTENING ADDENDA
It is evident that maximising β and therefore dilating the disturbance set is equivalent to
minimising δ. Hence, the policy max-dist is found by solving the convex optimisation problem
Pmax-dist(N) =

min
P ,{Z( · )},{Z¯( · )} δ
subject to:
(A.6), (A.7)
∃{Z(j) ≥ 0a×q}, {Z¯(j) ≥ 0a×r} :
s satisfies (5.17) and (5.18)
(A.8)
A.2 VH-ROA Plots for Different Policies
Presented below are plots for the full VH-ROA corresponding to the CT policies compared in
Section 5.5.
x1
x 2
-15 -10 -5 0 5 10 15
-5
-4
-3
-2
-1
0
1
2
3
4
5
N=4
N=3
N=2
N=1
Terminal set
Figure A.1: Nominal VH-ROA for T1
140
A.2. VH-ROA PLOTS FOR DIFFERENT POLICIES
x1
x 2
-8 -6 -4 -2 0 2 4 6 8
-5
-4
-3
-2
-1
0
1
2
3
4
5
N=4
N=3
N=2
N=1
Terminal set
Figure A.2: VH-ROA for T1 using max-vol
x1
x 2
-6 -4 -2 0 2 4 6
-5
-4
-3
-2
-1
0
1
2
3
4
5
N=4
N=3
N=2
N=1
Terminal set
Figure A.3: VH-ROA for T1 using max-dist
141
A. CONSTRAINT TIGHTENING ADDENDA
x1
x 2
-5 -4 -3 -2 -1 0 1 2 3 4 5
-4
-3
-2
-1
0
1
2
3
4
N=4
N=3
N=2
N=1
Terminal set
Figure A.4: VH-ROA for T1 using 3-step nil
x1
x 2
-4 -3 -2 -1 0 1 2 3 4
-4
-3
-2
-1
0
1
2
3
4 N=4
N=3
N=2
N=1
Terminal set
Figure A.5: VH-ROA for T1 using 2-step nil
142
A.2. VH-ROA PLOTS FOR DIFFERENT POLICIES
x1
x 2
-10 -5 0 5 10
-5
-4
-3
-2
-1
0
1
2
3
4
5
N=4
N=3
N=2
N=1
Terminal set
Figure A.6: Nominal VH-ROA for T2
x1
x 2
-5 0 5 10 15
-4
-3
-2
-1
0
1
2
3
4
5
N=4
N=3
N=2
N=1
Terminal set
Figure A.7: VH-ROA for T2 using max-vol
143
A. CONSTRAINT TIGHTENING ADDENDA
x1
x 2
-5 0 5 10 15
-2
-1
0
1
2
3
4
5
N=4
N=3
N=2
N=1
Terminal set
Figure A.8: VH-ROA for T2 using max-dist
x1
x 2
-5 0 5 10 15
-1
0
1
2
3
4
5
N=4
N=3
N=2
N=1
Terminal set
Figure A.9: VH-ROA for T2 using 3-step-nil
144
A.2. VH-ROA PLOTS FOR DIFFERENT POLICIES
x1
x 2
-5 0 5 10 15
-1
0
1
3
4
5
N=4
N=3
N=2
N=1
Terminal set
Figure A.10: VH-ROA for T2 using 2-step-nil
145
146
Appendix B
Shovel Model & Controller Details
B.1 Equations of Motion
This section contains the symbolic equations of motion in as derived in Subsection 6.2.2. Recall
that they can be expressed in the compact form
D(q, mp)q¨ + C(q, q˙, mp, m˙p) + G(q, mp) +Λq˙ = Φ. (B.1)
This equation can be written explicitly in terms of components as[
D11 D12
D21 D22
] [
r¨
θ¨
]
+
[
C1
C2
]
+
[
G1
G2
]
+
[
Λ11 Λ12
Λ21 Λ22
] [
r˙
θ˙
]
=
[
Fr
τθ
]
. (B.2)
In terms of the kinematic and dynamic parameters used in Section 6.2, the components of the
dynamic matrices are given by
D11 = M + mp (B.3)
D12 =
(
−r5 sin(1θ5)− 1ykcm
)
mp −M1y′cm (B.4)
D21 = D12 (B.5)
D22 =
(
r2 + 1y2kcm + r25 + 2r5 sin(1θ5)1ykcm + 2r1xkcm + 1x
2
kcm + 2r5 cos(
1θ5)r
+ 2r5 cos(1θ5)1xkcm
)
mp +
(
1y′2cm + r2 + 2r(1x′cm − Lc) + (1x′cm − Lc)2
)
M + J
(B.6)
C1 =
(
r˙− θ˙r5 sin(1θ5)− θ˙1ykcm
)
m˙p −
(
r5 cos(1θ5)θ˙2 + rθ˙2 + 1xkcmθ˙2
)
mp
−
(
rθ˙2 + (1x′cm − Lc)θ˙2
)
M
(B.7)
C2 =
(
θ˙1y2kcm − r˙1ykcm + θ˙r2 + θ˙r25 + 2θ˙r5 sin(1θ5)1ykcm − r˙r5 sin(1θ5) + 2θ˙r1xkcm
+ 2θ˙r5 cos(1θ5)1xkcm + 2θ˙r5 cos(1θ5)r + θ˙1x
2
kcm)
)
m˙p
+
(
2r˙θ˙1xkcm + 2r˙θ˙r5 cos(1θ5) + 2r˙θ˙r
)
mp +
(
2r˙θ˙(1x′cm − Lc) + 2r˙θ˙r
)
M
(B.8)
G1 = g
(
M sin(θ) + mp sin(θ)
)
(B.9)
147
B. SHOVEL MODEL & CONTROLLER DETAILS
G2 = g
((
r5 cos(θ + 1θ5)− sin(θ)1ykcm + cos(θ)1xkcm + cos(θ)r
)
mp
+
(
− sin(θ)1y′cm + cos(θ)r + cos(θ)(1x′cm − Lc)
)
M
)
.
(B.10)
The components of the frictional matrix Λ are given in (6.28).
B.2 Controller Parameters
The pre-stabilising controller gain is given by
K = 1× 107
[
0.5113 0.1593 0.7160 −0.1016 0.0034 −0.0004
0.0052 2.0962 0.0368 2.8854 −0.0001 0.0556
]
. (B.11)
The identified, downsampled linear system is given by the state space matrices
A =

−3.5326 1.6879 −5.5772 6.9162 0.0054 −0.0186 0.0271
−0.1732 0.5527 −0.2599 1.0306 0.0004 −0.0007 0.0019
6.5103 −5.5942 10.1230 −16.0435 −0.0103 0.0293 −0.0586
−1.8555 0.1872 −2.6759 3.9571 0.0030 −0.0079 0.0153
381.4447 −1942.8321 580.8633 −2955.5991 −0.5684 1.9160 −10.7236
−1170.9576 1168.1313 −1822.8751 3206.4945 1.9817 −4.6376 11.3804
488.6781 138.6933 700.0382 −773.7160 −0.7584 2.0613 −2.8468

(B.12)
B =

0.1562 0.0208
0.0007 0.5822
0.3149 0.0172
−0.0137 1.0105
45.1662 450.6227
−141.9427 −76.3676
−4.3556 −438.5990

(B.13)
C =
 I7[ −7.8467 15.5957 −12.5835 27.0824 0.0211 −0.0371 0.1016
−10.7536 20.2900 −14.8304 23.6959 0.0145 −0.0428 0.0864
] (B.14)
D = 09×2 (B.15)
148