Distributed MPC for PWA Systems Based on Switching ADMM
Samuel Mallick Azita Dabiri Bart De Schutter

	















Distributed MPC for PWA Systems Based on Switching ADMM
    Samuel Mallick, Azita Dabiri, and Bart De Schutter, Fellow, IEEE
	
	This paper is part of a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon
		2020 research and innovation programme (Grant agreement No. 101018826
		- CLariNet).
	All authors are affiliated with Delft Center for Systems and Control, Delft
		University of Technology, Delft, The Netherlands (e-mail: {s.h.mallick,a.dabiri,b.deschutter}@tudelft.nl).
    April 29, 2024
==============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





	This paper presents a novel approach for distributed model predictive control (MPC) for piecewise affine (PWA) systems.
	Existing approaches rely on solving mixed-integer optimization problems, requiring significant computation power or time.
	We propose a distributed MPC scheme that requires solving only convex optimization problems.
	The key contribution is a novel method, based on the alternating direction method of multipliers, for solving the non-convex optimal control problem that arises due to the PWA dynamics.
	We present a distributed MPC scheme, leveraging this method, that explicitly accounts for the coupling between subsystems by reaching agreement on the values of coupled states.
	Stability and recursive feasibility are shown under additional assumptions on the underlying system.
	Two numerical examples are provided, in which the proposed controller is shown to significantly improve the CPU time and closed-loop performance over existing state-of-the-art approaches.


	Distributed model predictive control, networked systems, piecewise affine systems, alternating direction method of multipliers (ADMM)



§ INTRODUCTION


Distributed control of large-scale hybrid systems is an important challenge in the field of systems and control <cit.>.
Hybrid systems are dynamical systems exhibiting both continuous and discrete dynamics.
They are powerful models for many real-world large-scale systems such as transportation networks <cit.>, power networks <cit.>, and water networks <cit.>.
In particular, piecewise affine (PWA) models are a popular class of hybrid system model due their suitability for analysis <cit.>, proven equivalence to other hybrid model classes <cit.>, and ability to approximate non-linear functions to arbitrary accuracy.

Distributed model predictive control (MPC) <cit.> is a promising candidate for control of large-scale PWA systems.
MPC is an optimization-based paradigm where control inputs are generated online by minimizing a sum of stage costs over a prediction horizon, subject to the system dynamics, and physical and operational constraints <cit.>.
Distributed MPC applies this idea to spatially separated or decomposed systems, where subsystem controllers solve, sometimes iteratively, local optimization problems, and where the coupling between connected subsystems is accounted for by communication <cit.>.
For linear systems, this approach can reconstruct the solution to a globally defined convex MPC controller through the use of iterative distributed optimization algorithms, e.g., the alternating direction method of multipliers (ADMM) <cit.>.
In this case, as the globally defined controller considers the global system, the effect of coupling between subsystems is accounted for and the stability and performance of the globally defined MPC controller carries through to the distributed implementation.
For PWA systems, however, the PWA dynamics result in a non-convex global optimization problem for which distributed optimization algorithms, such as ADMM, are not guaranteed to recover the global optimum, nor even to find a feasible solution.
As such, the effects of coupling between subsystems may cause instability or negatively affect the performance of the resulting distributed controller.
Furthermore, the local optimization problems are also non-convex, and the computational burden of solving them iteratively can become prohibitive, particularly as their computational complexity grows exponentially with the size of the prediction horizon <cit.>.
Therefore, two key challenges for distributed MPC for PWA systems are: how to reduce the online computational burden, and how to account for the effects of coupling between subsystems.

The few existing approaches for distributed MPC for PWA systems rely on converting the PWA dynamics into mixed logical dynamical (MLD) form, converting the MPC optimization problems into mixed-integer quadratic problems (MIQPs) <cit.>.
Distributed control is then conducted by the subsystems solving local MIQPs, iteratively in parallel <cit.>, or in a sequential order <cit.>, with communication of the local solutions between subsystems handling the coupling.
However, these approaches suffer from the computational complexity of solving multiple MIQPs.
To reduce the online computational burden, <cit.> proposes to select values for the integers in the MIQPs with heuristics, using convex distributed MPC techniques to then solve for the values of the continuous variables.
However, the heuristic choice of integer solutions is in general highly sub-optimal.
Perhaps the most theoretically complete approach was proposed recently in <cit.>, where distributed MPC for PWA systems is presented using robustly controllable sets.
The computational complexity of the local MIQPs is reduced by considering a prediction horizon of length one, and convergence to a terminal set is proven.
However, the coupling between subsystems is handled using robustness techniques, where subsystems consider coupling effects as disturbances, rather than agreeing on the values of the shared states.
This adds conservativeness, reducing performance and limiting the approach to systems with weak coupling.

In light of these challenges, this paper presents a novel approach for distributed MPC for PWA systems that systematically handles the coupling between subsystems and relies only on convex optimization.
The key contribution, enabling the approach, is a novel distributed method for solving the global non-convex optimization problem arising in the globally formulated MPC controller.
This method is based on ADMM, and leverages the underlying PWA dynamics of the system to solve the non-convex optimization problem distributively using only convex optimization.
Furthermore, the coupling between subsystems is handled explicitly, with the distributed solutions provably converging to points where the values of all shared states are agreed upon.
As such, the proposed method can improve upon the performance of existing controllers that heuristically handle coupling between subsystems, and presents a significantly lower computational burden than approaches based on mixed-integer optimization.

This paper is organized as follows.
Section <ref> gives the preliminaries, formalizing the systems considered and giving background on centralized and distributed MPC.
Section <ref> presents the switching ADMM procedure and the distributed MPC controller, with the stabilizing properties proved in Section <ref>.
In Section <ref> two numerical examples demonstrate the effectiveness of the approach, and, finally, Section <ref> concludes the paper.



§ PRELIMINARIES




 §.§ Notation

Define the following sets, all subsets of ℕ, as ℳ = {1,…,M}, ℒ_i = {1, …, L_i}, 𝒦 = {0,…,N-1}, and 𝒦' = {0,…,N}.
We use the counter t to represent closed-loop time steps and k for time steps within an MPC prediction horizon.
A vector that stacks the vectors x_i, i ∈ℳ, in one column vector is denoted col_i∈ℳ(x_i).
Define the closure of set 𝒫 as 𝒫.
For brevity we write (A)^⊤ B (A) as (⋆)^⊤ B (A).



 §.§ Distributed Piecewise Affine Systems

We consider a global system composed of M subsystems where each subsystem i ∈ℳ has a state x_i ∈ℝ^n_i and an input u_i ∈ℝ^m_i. 
Both states and inputs are subject to convex constraints

    x_i ∈𝒳_i,    u_i ∈𝒰_i,  ∀ i ∈ℳ.

The graph 𝒢 = (ℳ, ℰ) defines a coupling topology where the edges ℰ are ordered pairs (i, j) indicating that subsystem i effects subsystem j through the dynamics, stage cost, constraints, or some combination thereof. 
Define the neighborhood of subsystem i as 𝒩_i = {j ∈M | (j, i) ∈ℰ, i ≠ j}.
A subsystem is not in its own neighborhood, i.e., (i,i) ∉ℰ.


	Subsystem i and subsystem j can communicate in a bidirectional manner if i ∈𝒩_j or j ∈𝒩_i.


The subsystems are governed by discrete-time piecewise affine (PWA) dynamics, with the state update equations f_i affine on a finite number of convex polytopes, {P_i^(l)}_l ∈ℒ_i, referred to as regions:

    x_i^+    = f_i(x_i, u_i, {x_j}_j ∈𝒩_i) 
       = A_i^(l) x_i + B_i^(l) u_i + c_i^(l)
          + ∑_j ∈𝒩_i A_ij^(l) x_j,   x_i ∈ P_i^(l), ∀ l ∈ℒ_i.

blackThe regions have non-overlapping interiors, i.e., int(P_i^(l)) ∩int(P_i^(l^')) = ∅ for l,l' ∈ℒ_i, l ≠ l', and form a partition of the state space, i.e., ⋃_l ∈ℒ_i P_i^(l) = 𝒳_i.
Define the global state and input, the aggregation of states and inputs for all subsystems, as x = col_i ∈ℳ(x_i) and u = col_i ∈ℳ(u_i).
Furthermore, define variables that gather states and inputs over the prediction horizon: x_i = (x_i^⊤(0),...,x_i^⊤(N))^⊤, u_i = (u_i^⊤(0),...,u_i^⊤(N-1))^⊤, x = (x(0)^⊤,...,x(N)^⊤)^⊤, u = (u(0)^⊤,...,u(N-1)^⊤)^⊤.



 §.§ Centralized MPC


Consider the global MPC controller defined by the following finite-horizon optimal control problem parametrized by x:

	
    𝒫(x):  min_x,u ∑_i ∈ℳ F_i( x_i, u_i, {x_j}_j ∈𝒩_i)
    		
    s.t.     ∀  i ∈ℳ: 
       x_i(k+1) = f_i(x_i(k), u_i(k), {x_j(k)}_j ∈𝒩_i),  ∀ k ∈𝒦
       h_i(x_i(k), {x_j(k)}_j ∈𝒩_i) ≤ 0,  ∀ k ∈𝒦'
       (x_i(k), u_i(k)) ∈𝒳_i ×𝒰_i,  ∀ k ∈𝒦' 
       x_i(0) = x_i,

	
with

    F_i   ( x_i, u_i, {x_j}_j ∈𝒩_i) = 
       ∑_k ∈𝒦ℓ_i(x_i(k), u_i(k), {x_j(k)}_j ∈𝒩_i) + V_f, i(x_i(N)).

The convex functions ℓ_i and V_f, i are stage and terminal costs, and the linear functions h_i define coupled convex inequality constraints.
The first entries in the optimal input sequences define feedback control laws u_MPC, i(x) = u^⋆_i(0).

Define the set of states x for which 𝒫 is feasible as 𝒳_0, and define a set of local control sequences {u_i}_i ∈ℳ to be globally feasible for x if the corresponding global control sequence u is feasible for 𝒫(x).
For a global state x ∈𝒳_0 and feasible control sequence u define the value as

    V(x, u) = ∑_i ∈ℳ F_i( x_i, u_i, {x_j}_j ∈𝒩_i),

where the state sequences x_i are generated by applying u to the dynamics (<ref>) starting from x.

The MPC scheme defined by 𝒫 considers general coupling between subsystems, i.e., in the costs, dynamics, and constraints.
Solving 𝒫 directly is undesirable as the complexity of the optimization problem grows with the number of subsystems.
Additionally, a centralized solution requires centralized coordination, which cannot preserve privacy of the local functions of subsystems, and presents a single point of failure.
It is therefore desirable to solve 𝒫 distributively, where each subsystem computes a local control input using only local information and communication with its neighbors.

	The approach presented in this paper extends easily to the case where the cost functions F_i also switch depending on the subsystem's PWA region.
For clarity of exposition we do not explicitly include this case in 𝒫, as the notation would then become unwieldy.




 §.§ Distributed MPC

By introducing copies of coupled states the problem 𝒫 can be reformulated for distributed optimization <cit.>.
Define an augmented state trajectory for each subsystem that includes copies of the state trajectories of neighboring subsystems

    x_i = (x_i^⊤, col^⊤_j ∈𝒩_i(x_j^(i)))^⊤,

where x_j^(i) is agent i's local copy of agent j's state trajectory.
Furthermore, define a global augmented state as x = col_i ∈ℳ(x_i).
An equivalent problem to 𝒫 is then

	
    𝒫_d(x):  min_x, u ∑_i ∈ℳ F_i( x_i, u_i, {x_j^(i)}_j ∈𝒩_i)
    		
    s.t.     ∀ i ∈ℳ: 
       x_i(k+1) = f_i(x_i(k), u_i(k), {x_j^(i)(k)}_j ∈𝒩_i),  ∀ k ∈𝒦
       h_i(x_i(k), {x_j^(i)(k)}_j ∈𝒩_i) ≤ 0,  ∀ k ∈𝒦' 
       (<ref>), (<ref>)
       x_j^(i)(k) = x_j(k),   j ∈𝒩_i,  ∀ k ∈𝒦'.


The equality constraint (<ref>) gives consistency between the true state trajectories and their copies and therefore equivalence between 𝒫 and 𝒫_d. 
In 𝒫_d the stage costs, dynamics, and inequality constraints can be separated across the local variables x_i and u_i.
The alternating direction method of multipliers (ADMM) <cit.> can then be applied to solve 𝒫_d distributively by dualizing the consistency constraint (<ref>), and having systems iteratively solve and communicate the solutions to local sub-problems derived from 𝒫_d.
However, due to the PWA dynamics (<ref>) 𝒫_d is non-convex, and ADMM is hence not guaranteed to converge to an optimal, or even feasible, solution.
Additionally, systems must solve non-convex local optimization problems iteratively, requiring significant computation time or power.
To address this, in the following section we introduce distributed MPC based on a switching ADMM procedure that, leveraging the structure of the PWA dynamics, can guarantee convergence to a feasible solution, and requires solving only convex optimization problems.



§ SWITCHING ADMM-BASED DISTRIBUTED MPC

In this section we introduce the proposed distributed MPC algorithm.
First, we show that the structure of 𝒫_d is piecewise-convex over a finite collection of polytopes.
This structure is then leveraged to formulate the switching ADMM procedure upon which the distributed MPC algorithm is based.



 §.§ Structure of 𝒫_d


The following analysis is inspired by the reverse transformation presented in <cit.>, where MPC for a single PWA system was considered.
Define concatenated decision variables as y_i = (x_i^⊤, u_i^⊤)^⊤ and y = (x^⊤, u^⊤)^⊤.
The feasible set for 𝒫_d(x) is then

    𝒴(x) = {y = (x^⊤, u^⊤)^⊤| (<ref>)­-(<ref>), ∀ i ∈ℳ}.

We now introduce the notion of switching sequences that specify a sequence of PWA regions over the prediction horizon.
Consider the switching sequence s_i ∈𝒮_i = {1,...,L_i}^N+1 that specifies the PWA regions for subsystem i as:

    s_i = (s_i(0), s_i(1), ..., s_i(N)) 
        x_i(k) ∈ P_i^(s_i(k)),  ∀ k ∈𝒦'.

and the corresponding time-varying linear dynamics

    x_i(k+1)    = f_i^(s_i)(x_i(k), u_i(k), {x_j^(i)(k)}_j ∈𝒩_i) 
       = A_i^(s_i(k)) x_i(k) + B_i^(s_i(k)) u_i(k) 
          + c_i^(s_i(k)) + ∑_j ∈𝒩_i A_ij^(s_i(k)) x_j^(i)(k),  ∀ k ∈𝒦.

A particular x_i is said to generate a switching sequence s_i if 

    x_i(k) ∈P_i^(s_i(k)),  ∀ k ∈𝒦'.

Note that more than one s_i can be generated by x_i if x_i(k) lies on the boundary of at least two PWA regions for some k, i.e., if x_i(k) ∈P_i^(a)∩P_i^(b) then two sequences are generated, one containing s_i(k) = a and the other s_i(k) = b.
For a given x_i, u_i, and {x_j^(i)}_j ∈𝒩_i subsystem i can evaluate the generated switching sequences by rolling out and branching its dynamics as in Algorithm <ref>.


Let s = (s_1, …, s_M) ∈𝒮, 𝒮 = 𝒮_1 ×…×𝒮_M, define a global switching sequence specifying local switching sequences for each subsystem.
This global s defines a new optimal control problem

	
    𝒫_s(x):     min_y ∑_i ∈ℳ F_i( x_i, u_i, {x_j^(i)}_j ∈𝒩_i)
    		
    s.t.     ∀  i ∈ℳ: 
       (<ref>)­-(<ref>), (<ref>)
       x_i(k) ∈ P_i^(s_i(k)),  ∀ k ∈𝒦',


where the time-varying dynamics (<ref>) replace the PWA dynamics and the PWA regions selected by s are enforced with the additional constraints (<ref>).
The feasible set for 𝒫_s(x) is

    𝒴_s(x) = {y = (x^⊤, u^⊤)^⊤| (<ref>) - (<ref>), (<ref>), (<ref>), ∀ i ∈ℳ},

and it is convex as the constraints (<ref>)­-(<ref>) are convex, the dynamics (<ref>) are linear, and the PWA regions in (<ref>) are convex.
Note that 𝒴_s is empty for many choices of s.
We now show that 𝒴 is the union of the sets 𝒴_s.

	Suppose x and s are given such that 𝒴_s(x) ≠∅. 
	Then for all y = (x^⊤, u^⊤)^⊤∈𝒴_s(x):
	
    f_i(x_i(k), u_i(k), {x_j^(i)(k)}_j ∈𝒩_i) 
          = f_i, s(x_i(k), u_i(k), {x_j^(i)(k)}_j ∈𝒩_i),
	
	∀ k ∈𝒦,∀ i ∈ℳ.


	Let x and s be given, and let y be an arbitrary element of 𝒴_s(x). 
	From (<ref>), x_i(k) ∈ P_i^(s_i(k)), ∀ k ∈𝒦, ∀ i ∈ℳ, so that, from (<ref>),
	
    f_i(x_i(k), u_i(k), {x_j^(i)(k)}_j ∈𝒩_i) 
       = (A_i^(l) x_i + B_i^(l) u_i + ∑_j ∈𝒩_i A_ij^(l))|_l = s_i(k)
        = f_i, s(x_i(k), u_i(k), {x_j^(i)(k)}_j ∈𝒩_i),

	for all k ∈𝒦 and i ∈ℳ.


	The feasible set for problem 𝒫_d(x) is the union of a finite number of convex sets, that are the feasible sets for the problems 𝒫_s(x), i.e.,
	
    𝒴(x) = ⋃_s ∈𝒮𝒴_s(x).



	i) Suppose y∈𝒴(x).
	Let s be a switching sequence generated by y, so that x_i(k) ∈ P_i^(s_i(k)), ∀ k ∈𝒦, ∀ i ∈ℳ.
	Hence, y satisfies (<ref>) in 𝒫_s(x).
	Furthermore, the dynamics in (<ref>) reduce to (<ref>).
	Finally, y satisfies (<ref>)­-(<ref>) by definition.
	Thus, y∈𝒴_s(x) and 𝒴(x) ⊆𝒴_s(x) ⊆⋃_s ∈𝒮𝒴_s(x).
	ii) Now suppose y∈⋃_s ∈𝒮𝒴_s(x).
	Then there exists s ∈𝒮^M such that y∈𝒴_s(x).
	By Lemma <ref>, for all y∈𝒴_s(x), (<ref>) is equivalent to the dynamics in (<ref>), and y satisfies (<ref>).
	Additionally, y satisfies (<ref>)­-(<ref>) by definition.
	Thus,  y∈𝒴(x) so that ⋃_s ∈𝒮𝒴_s(x) ⊆𝒴(x).
	Equation (<ref>) follows.


Proposition <ref> reveals the structure of 𝒫_d. 
It is piecewise convex, with each convex region associated with a global switching sequence s and being defined over the set 𝒴_s.
It follows that, in theory, 𝒫_d can be solved by solving 𝒫_s for each s ∈𝒮, and taking the solution that gives the lowest cost.
However, in practice, this is not feasible as |S| grows exponentially with M and N.



 §.§ Switching ADMM Procedure

Rather than solving every convex piece 𝒫_s, we propose to consider one piece at a time.
Applying ADMM to one 𝒫_s will converge to the optimum of that convex problem.
However, during the iterations it can be beneficial to switch to a new 𝒫_s that gives a lower cost.
The proposed idea is then to allow subsystems, during the ADMM iterations, to locally change switching sequences s_i, changing the global s and 𝒫_s, thereby moving from one convex piece to another convex piece of 𝒫_d, until a local minimum is found.
In this section we formalize this idea.

First, we make the ADMM procedure for solving 𝒫_s explicit.
Define the local constraint set containing the constraints independent of s_i as

    𝒴_i(x_i) =    {y_i | (<ref>), (<ref>)),

and the local constraint set containing the constraints dependent on s_i as

    𝒴_i^(s_i) =    {y_i |(<ref>), (<ref>)}.

Introduce additional copies of each state trajectory z = (z_1^⊤, z_2^⊤, ...., z_M^⊤)^⊤, with the relevant components of z for subsystem i denoted z_i = (z_i^⊤, col^⊤_j ∈𝒩_i(z_j))^⊤.
The problem 𝒫_s(x) can then be reformulated as

    min_{y_i}_i ∈ℳ∑_i ∈ℳ   F_i( x_i, u_i, {x_j^(i)}_j ∈𝒩_i) 
    s.t.     ∀  i ∈ℳ:
         y_i ∈𝒴_i(x_i) ∪𝒴_i^(s_i)
       x_i = z_i,

with the final constraint forcing agreement across all coupled states.
Recall that in (<ref>) the switching sequences s_i are not optimization variables and have been prespecified.

ADMM solves (<ref>) with the following distributed iterations <cit.>:

	
	
    y_i^τ+1   = _y_i ∈𝒴_i ∪𝒴_i^(s_i)   F_i( x_i, u_i, {x_j^(i)}_j ∈𝒩_i) + (λ_i^τ)^⊤x_i 
             + ρ/2x_i - z_i^τ_2^2 
    z_i^τ+1   = 1/|𝒩_i|+1 (x_i^τ+1 + ∑_j ∈𝒩_ix_i^(j), τ+1) 
    λ_i^τ+1   = λ_i^τ + ρ(x_i^τ+1 - z_i^τ+1),


with λ_i the Lagrange multiplier for subsystem i and ρ > 0 a penalty parameter.
This is a distributed message passing algorithm where, in step (<ref>) each subsystem solves a local optimization problem, in step (<ref>) the local solutions are communicated and averaged between neighboring subsystems, and in step (<ref>) the subsystems update their multipliers locally.
Define the residual r(y) as the summed disagreement on the values of shared states

    r(y) = ∑_i ∈ℳ∑_j ∈𝒩_ix_j^(i) - x_j_2.
 
As τ→∞ the local solutions converge to the minimizers of 𝒫_s, i.e., y_i^τ→y_i^∗, ∀ i ∈ℳ, and the shared states converge to agreement <cit.>

    x_j^(i), τ - x_j^τ→ 0,  ∀ j ∈𝒩_i,  ∀ i ∈ℳ,

i.e., r(y^τ) → 0.
Additionally, as the optimizer of 𝒫_s is trivially a feasible point of 𝒴_s, and 𝒴_s ⊂𝒴, we have that as τ→∞

    y^τ∈𝒴.

In the following we will refer to any point where (<ref>) and (<ref>) hold as a consensus point.

Our key insight is that new local switching sequences can be identified when the local solutions to (<ref>) are pushed to the boundary of the feasible set 𝒴_i^(s_i).
Indeed, in the case where the current 𝒫_s does not contain a local minimum of 𝒫_d this is the expected behavior, as the minimum of 𝒫_s is then on the boundary of its feasible set.
When this occurs, subsystems change their local switching sequences s_i, which globally changes s and 𝒫_s.
In this way, the non-convex problem 𝒫_d is explored by jumping between the different approximations 𝒫_s, with the jumping determined distributively by intermediate local solutions to (<ref>).
Once local solutions are no longer pushed to boundaries of 𝒴_i^(s_i) the system-wide 𝒫_s no longer changes, and the ADMM iterations will converge, as 𝒫_s is convex, to a local minimum of 𝒫_d and a consensus point.

Let us elaborate with the representative visual example in Fig.  <ref> that depicts two systems, each with two PWA regions P_i^(1) and P_i^(2) and current states x_1 ∈ P_1^(1) and x_2 ∈ P_2^(1).
To facilitate visualization consider scalar states and inputs u_i, x_i ∈ℝ and a horizon of N = 1, such that the control input decision variables for each subsystem are scalar and the switching sequences are of length two. 
There are 16 possible global switching sequences 

    𝒮 = {s = (s_1, s_2) | s_1, s_2 ∈{(1, 1), (1, 2), (2, 1), (2, 2)}}.

However, as x_1 ∈ P_1^(1) and x_2 ∈ P_2^(1), only switching sequences with s_1(0) = s_2(0) = 1 correspond to a non-empty 𝒫_s.
The top plot of Fig. <ref> shows, over the control input decision space, the four non-empty convex pieces, {𝒫_s}_s ∈𝒮̂, 𝒮̂ = {(s_1, s_2) | s_1, s_2 ∈{(1, 1), (1, 2)}}, constituting 𝒫_d.
The bottom plot shows the local decision variables x_1(1) and x_2(1) and the PWA partitions.   
An initial guess for u_1(0) and u_2(0) generates x_1(1) ∈ P_1^(1) and x_2(1) ∈ P_2^(1) (blue dot), and thus s_1 = s_2 = (1, 1).
Each subsystem sets the constraint set 𝒴_i^(s_i) that defines the problem 𝒫_s (yellow region).
The subsystems commence the ADMM iterations performing one round of the steps in (<ref>).
The output of the local minimization yields new solutions (orange dot), and subsystem 2 identifies that this local solution is on the boundary of P_2^(1) and P_2^(2) using Algorithm <ref>. 
Subsystem 2 then changes its switching sequence to s_2 = (1, 2), changing the local constraint set 𝒴_2, s_2, and changing the global convex problem to a new 𝒫_s (purple region).
The ADMM iterations continue, now with subsystem 2 having changed its local minimization problem (<ref>) by changing 𝒴_2^(s_2). 
At the fourth iteration (purple dot), subsystem 1 identifies a new switching sequence s_1 = (1, 2), and changes 𝒴_1^(s_1), changing the global problem to the red region. 
Finally, as the local solutions are no longer pushed to the boundaries of the local constraint sets 𝒴_i^(s_i), the iterations converge to the minimum of the convex piece (red region), a local minimum for 𝒫_d.


Generating local switching sequences with Algorithm <ref> requires a set of local copies {x_j^(i)}_j ∈𝒩_i.
While these are an output of the local minimization step (<ref>), for the initial switching sequence, i.e., determining the yellow region from the blue dot in Fig.  <ref>, these values can be attained with the distributed procedure outlined in Algorithm <ref> involving N steps of rolling out local dynamics and neighbor-to-neighbor communications.


Algorithm <ref> formalizes the switching ADMM procedure.
Note that each subsystem requires a local control sequence guess to generate the initial local switching sequences.
These guesses are discussed in Sections <ref> and <ref>.

We introduce a cut-off number of iterations T_cut, after which the switching sequences s_i no longer change.
This protects against the edge case where local solutions pull each other into new switching sequences with equal force, resulting in continuous oscillation between two convex regions of 𝒫_d.
We note that in our experiments this case is observed only rarely, with the iterations otherwise converging to a problem 𝒫_s that contains a local minimum of 𝒫_d.

We now give the convergence result for Algorithm <ref>.
Following the introduction of T_cut an additional assumption is required to prove convergence of the algorithm.

	
	For all x ∈𝒳_0 there does not exist a global switching sequence s = (s_1, …, s_M) such that 𝒴_i(x_i) ∪𝒴_i^(s_i)≠∅,  ∀ i ∈ℳ and 𝒴_s(x) = ∅.

Assumption <ref> protects against the edge case where, when τ = T_cut, the local systems are in a configuration of local switching sequences that represent a 𝒫_s with an empty feasible set.
 
	For all x ∈𝒳_0, under Assumption <ref>, as T_→∞, the outputs of Algorithm <ref>, {y_i^T_}_i ∈ℳ, converge to a consensus point, i.e., r(y^T_) = 0 and y^T_∈𝒴.


	For iterates τ≥ T_cut the local sets 𝒴_i^(s_i) are fixed and, by Assumption <ref>, the corresponding 𝒫_s has a non-empty feasible set. Therefore, as 𝒫_s is convex, the standard ADMM iterations (<ref>) converge to a consensus point <cit.>.


	Assumption <ref> is an assumption on the strength of the coupling between subsystems.
	Providing conditions to easily verify Assumption <ref> is beyond the scope of the current paper.
	An alternative approach, avoiding Assumption <ref>, would be to permit subsystems to switch their local switching sequence only after the new global convex region has been checked for feasibility.
	This is left to future work.




 §.§ Distributed MPC Algorithm
 
Leveraging the switching ADMM procedure we now present the distributed MPC algorithm in Algorithm <ref>.
In step <ref> an initial globally feasible set of local guesses is generated.
In Section <ref> it is shown how subsystems can locally generate these guesses using a terminal control law.
Here, however, we highlight that these initial guesses can often be generated using heuristics specific to the application.
Additionally, a common heuristic is to use the shifted solution from the previous time step, i.e., u_i = (u̅_i^⊤(1),…,u̅_i^⊤(N-2), u̅_i^⊤(N-2))^⊤, where u̅_i = (u̅_i^⊤(0),…,u̅_i^⊤(N-1))^⊤ is the solution of the previous time step.
Finally, state constraints can be made soft constraints through the introduction of slack variables, penalized in the objective function, to ensure any initial feasible guess is globally feasible.

As Algorithm <ref> generates, in general, a local minimum of 𝒫_d, the quality of the solution could be improved by running it in parallel with different initial guesses, introducing an effective multi-start approach.
This is of course limited by the communication and computation limits of the subsystems.


	While the convergence of ADMM, and hence Algorithm <ref>, is asymptotic, in practice ADMM often converges to acceptable accuracy within few iterations <cit.>, motivating the use of a finite T_ADMM.
	Furthermore, if errors induced by finite termination were to become an issue, existing work that ensures feasibility and stability under finitely terminated distributed optimization <cit.> can be added to our approach.
	 




§ STABILITY AND RECURSIVE FEASIBILITY

In this section we give a stable and recursively feasible extension to Algorithm <ref> under additional assumptions on the system.
Given the local properties of Algorithm <ref>, we leverage a dual-mode approach to stability <cit.>, where terminal constraints are added to the MPC scheme and subsystems use stable switching linear controllers when they enter a terminal set.
This introduces a complex trade-off between the size of the terminal set, the length of the prediction horizon, and closed-loop performance.
For a larger terminal set the systems use linear control laws for a larger portion of the state-space, using the optimization-based MPC control law less and foreseeably sacrificing performance.
With a smaller terminal set, on the other hand, the MPC controller is used more, but requires a longer horizon to satisfy the terminal constraint.
We note that, as Algorithm <ref> involves only convex programs, long horizons present less computational demand than for approaches that rely on mixed-integer programming.



 §.§ Assumptions

Define ℒ̂_i = {l ∈ℒ_i | 0 ∈P_i^(l)} the set of subsystems i's PWA regions containing or touching the origin.
We introduce the following extra assumptions on the systems considered.
[<cit.>] 
	For all i ∈ℳ the matrix pairs (A_i^(l), B_i^(l)),  ∀ l ∈ℒ_i, are controllable.
	The origin is an equilibrium point with u_i = 0 and x_j = 0,  ∀ j ∈𝒩_i, i.e., c_i^(l) = 0,  ∀ l ∈ℒ̂_i.
	There exist for all l ∈ℒ̂_i a gain K_i^(l) such that A_i^(l) + B_i^(l)K_i^(l) is a Schur matrix. 

Assumption <ref> says that, for each subsystem, the origin is an equilibrium of the dynamics in all PWA regions that touch it, and that there exists a stabilizing linear feedback law within each of these regions.
 
	For all i ∈ℳ there exists a non-empty terminal set 𝒳_T, i⊆⋃_l ∈ℒ̂_i P_i^(l), containing the origin, such that ∀ x_i ∈𝒳_T, i we have f_i(x_i, u_i, {x_j}_j ∈𝒩_i) ∈𝒳_T, i,  ∀ x_j ∈𝒳_T, j,  ∀ j ∈𝒩_i, with u_i = K_i^(l) x_i and K_i^(l) x_i ∈𝒰_i. 

Define the global terminal set to be 𝒳_T = {x | x_i ∈𝒳_T, i, ∀ i ∈ℳ}.
Assumption <ref> guarantees that there exist local terminal sets that are positively invariant under the switching terminal controllers and that are robust to the effects of coupling when neighboring subsystems are also within their terminal sets.
This is a less restrictive assumption than that in <cit.>, where the terminal sets are assumed to be robust to the effects of coupling for the entire state space of neighboring systems 𝒳_j.
In contrast to <cit.>, this less restrictive assumption is possible due to the proposed approach explicitly giving agreement on the values of shared states over the prediction horizon, rendering the approach applicable to systems with stronger coupling.
Note that <cit.> proposes a method to compute the terminal sets required for the controller in <cit.>.
These sets trivially satisfy Assumption <ref> as 𝒳_T, j⊂𝒳_j.
blackIn Appendix <ref> we present a modification to the procedure given in <cit.> to calculate less restrictive sets satisfying Assumption <ref>.
Additionally, we provide a formal proof, missing in <cit.>, for the forward invariance property of the resulting sets.
Clearly, the sets required by <cit.> do not exist for systems with unbounded state spaces, while the sets in Assumption <ref> can still exist.
 
	The stage and terminal costs are of the form ℓ_i(x_i, u_i, {x_j}_j ∈𝒩_i) = x_i^⊤ Q_i x_i + u_i^⊤ R u_i + ∑_j ∈𝒩_i x_j^⊤ Q_ij x_j and V_f, i(x) = x_i^⊤Φ_i x_i for positive definite Q_i, R_i, and Φ_i. 
	The terminal costs V_f, i and controllers K_i^(l) are such that for all x ∈𝒳_T,
	
    (∑_i ∈ℳ   V_f, i( f_i(x_i, u_i, {x_j}_j ∈𝒩_i) ) - V_f, i(x_i) 
            + ℓ_i(x_i, u_i, {x_j}_j ∈𝒩_i) ) ≤ 0,

	where u_i = K_i^(l) x_i, x_i ∈ P_i^(l).

Assumption <ref> says that the terminal cost matrices Φ_i result in ∑_i ∈ℳV_f, i being a global common quadratic Lyapunov function for all possible combinations of local PWA regions within the terminal set.
The computation of such terminal costs has been covered for single PWA systems in <cit.> and for distributed linear systems in <cit.>.
In Appendix <ref> we outline an LMI approach to constructing terminal costs V_f, i that satisfy Assumption <ref>.
blackWhile Assumption <ref> may be restrictive if the origin is on the boundary of many PWA regions, note that we could restrict the local terminal sets to a single PWA region, 𝒳_T, i⊆ P_i^(l_i), l_i ∈ℒ̂_i, such that finding a terminal cost satisfying Assumption <ref> reduces to finding a Lyapunov function for the linear system defined by the PWA regions {P_i^(l_i)}_i ∈ℳ, rather than a common Lyapnuov function over arbitrary switching.
Again we note that, as Algorithm <ref> involves only convex programs, long horizons are computationally tractable and smaller terminal sets are therefore viable.

	An initial globally feasible set of control sequences is assumed to be known for the first time step.

Assumption <ref> is common in distributed MPC <cit.>.



 §.§ Stabilizing Distributed MPC Algorithm

With a slight abuse of notation, we redefine the optimal control problems 𝒫_d(x) and 𝒫_s(x) by modifying the feasible sets 𝒴 in (<ref>) and 𝒴_i in (<ref>) to include terminal constraints

    𝒴(x)    = {y| (<ref>) - (<ref>), x_i(N) ∈𝒳_T, i, ∀ i ∈ℳ}
    𝒴_i(x_i)    = {y_i | (<ref>), (<ref>), x_i(N) ∈𝒳_T, i},

naturally redefining 𝒳_0 = {x   |  𝒴(x) ≠∅}.
We then introduce the stabilizing distributed MPC algorithm in Algorithm <ref>.

When all subsystems are in the terminal regions the switching linear feedback control laws are used.
Otherwise, for Algorithm <ref>, subsystems have local access to a feasible initial guess through shifting the previous control solution and adding the terminal controller as the final element of the initial feasible guess.

We now prove the recursive feasibility and the closed-loop stability of Algorithm <ref>.
 
	For an initial state x ∈𝒳_0 ∖𝒳_T, at time step t = 0, a globally feasible set of initial control guesses {u_i}_i ∈ℳ, and generating control inputs via Algorithm <ref>, 𝒫_d(x) is feasible for all time steps t ≥ 0 in which x ∉𝒳_T.


	The proof works by induction.
	Assume that, for state x ∉𝒳_T, at time step t, 𝒫_d(x_t) is feasible, and that subsystems have globally feasible initial control guesses u_i, ∀ i ∈ℳ.
	By Proposition <ref>, with feasible initial guesses, the outputs y_i = (x_i^⊤, u_i^⊤)^⊤ of Algorithm <ref> are feasible for 𝒫_d(x), and therefore x_i(N) ∈𝒳_T, i, ∀ i ∈ℳ.
	The subsystems apply local inputs u_i(0), propagating the state dynamics as x_i^+ = f_i(x_i, u_i(0), {x_j}_j ∈𝒩_i) and, as by Proposition <ref> x_j^(i) = x_j, ∀ j ∈𝒩_i, the new states coincide with the first predicted state, i.e., x_i^+ = x_i(1).
	At the next time step t+1 at state x^+ the initial guesses are  u^+_i = (u_i^⊤(1),…,u_i^⊤(N-1), (K_i^(l)x_i(N))^⊤)^⊤, x_i(N) ∈ P_i^(l)∩𝒳_T, i,   i ∈ℳ, and are globally feasible as 𝒳_T, i is forward invariant and robust to coupling under the terminal controllers when x_j ∈𝒳_T, j, ∀ j ∈𝒩_i.
	Therefore 𝒫_d(x^+) is feasible at time step t+1.
	
	Now, for time step t = 0 and initial state x ∈𝒳_0 ∖𝒳_T, 𝒫_d(x) is feasible as x ∈𝒳_0, and the initial guesses {u_i}_i ∈ℳ for Algorithm <ref> are globally feasible.
	It follows then by induction that 𝒫_d(x) is feasible for all time steps t > 0 in which x ∉𝒳_T.


	Theorem <ref> considers only states x ∉𝒳_T as 𝒫_d is not relevant in Algorithm <ref> when x ∈𝒳_T.


	The control law generated by Algorithm <ref> is exponentially stabilizing with a region of attraction 𝒳_0.


	For each subsystem i take u_i as the output from Algorithm <ref> at time step t with global state x, û_i = (u_i^⊤(1),…,u_i^⊤(N-1), (K_i^(l)x_i(N))^⊤)^⊤, x_i(N) ∈ P_i^(l) the initial guess at the next time step t+1, constructed from u_i, x^+ the global state at time step t+1 after applying local inputs u_i(0), and u^+_i as the output from Algorithm <ref> at time step t+1.
	Taking u, û, and u^+ naturally as the globally concatenated variables, we have from Assumption <ref>, by standard arguments,
	
    V(x^+, u^+) - V(x, u) ≤ V(x^+, û^+) - V(x, u) 
         ≤∑_i ∈ℳℓ_i (x_i, u(0), {x_j}_j ∈𝒩_i)

	for all x ∈𝒳_0 ∖𝒳_T, with V defined in (<ref>).
	Assumption <ref> on ℓ_i guarantee the existence of a class-𝒦 function β(·) such that 
	
    ∑_i ∈ℳℓ_i(x_i, u_i, {x_j}_j ∈𝒩_i) ≥β((x, u)),

	for all x ∉𝒳_T, and for all u ∈𝒰_1 ×…×𝒰_M.
	As 𝒳_T contains the origin there exists b > 0 such that V(x^+, u^+) - V(x, u) ≤ b for all x ∉𝒳_T.
	There therefore exists, for any initial state x ∈𝒳_0, a finite time t̅ at which the global state enters the global terminal set x ∈𝒳_T <cit.>. 
	As all local systems use the terminal controllers u_i = K_i^(l)x_i, x_i ∈ P_i^(l) for x_i ∈𝒳_T, i, and the local terminal sets are robustly positively invariant under this terminal control, we have x ∈𝒳_T for all t ≥t̅.
	Finally, by Assumption <ref>, under the switching linear controllers there exists a global common quadratic Lyapunov function for all possible combinations of local PWA regions when the global state x ∈𝒳_T.
	A common quadratic Lyapunov function implies exponential stability under arbitrary switching <cit.>, hence the global system is exponentially stable under the terminal controllers.
	Therefore as, for any initial state x ∈𝒳_0, the global state enters the terminal region 𝒳_T in finite time and is exponentially stable within 𝒳_T, the origin is exponentially stable under Algorithm <ref> with regions of attraction 𝒳_0.




§ ILLUSTRATIVE EXAMPLES


In this section we provide two examples demonstrating our approach.
The first considers the network from <cit.>, satisfying Assumptions <ref>­-<ref>.
This example demonstrates the stabilizing properties of Algorithm <ref> and compares the approach to that of <cit.>.
The second example considers the more realistic control challenge of a platoon of vehicles, proposed as a distributed hybrid MPC benchmark in <cit.>, and demonstrates the efficacy of Algorithm <ref> as a practical control scheme.
All examples are simulated on an 11th Gen Intel laptop with four i7 cores, 3.00GHz clock speed, and 16Gb of RAM.
Source code for the simulations can be found at <https://github.com/SamuelMallick/distributed-mpc-pwa, stable-dmpc-pwa/tree/paper_2024, hybrid-vehicle-platoon/tree/paper-2024>.
 


 §.§ Stabilizing Example

In the following we refer to Algorithm <ref> as SwA and the approach based on robust controllable sets presented in <cit.> as RCS.
Consider the network of three PWA systems, from <cit.>, with identical dynamics defined over L = 4 regions, as shown in Fig.  <ref>.
The systems' dynamics parameters, for i = 1,2,3, are: 

    A_i^(1)   = A_i^(3) = [ 0.6324 0.2785; 0.0975 0.5469 ]
       A_i^(2)   = A_i^(4) = [ 0.6555 0.7060; 0.1712 0.0318 ]
    
    		B_i^(1)   = B_i^(2) = B_i^(3) = B_i^(4) = [ 1; 0 ],

with state and action constraints 𝒳_i = {x_i | |x_i|_∞≤ 20} and 𝒰_i = {u_i | |u_i| ≤ 3}, and stage costs ℓ_i = x_i^⊤ Q x_i + u_i^⊤ R u_i with

    Q_i = [ 2 0; 0 2 ],     R = 0.2.
 
The coupling in the network is 

    A_12 = A_21 = A_23 = A_31 = 2 · 10^-3·[ 1 0; 0 1 ],

with all other coupling matrices zero.
For a fair comparison[In <cit.> a small disturbance is added to the dynamics and is accounted for in the computation of the robustly controllable sets. The approach of <cit.> requires no extra mechanism to account for the disturbance as the coupling between subsystems is also considered as a local disturbance. For a fair comparison we recompute all the sets required for the RCS approach with no additional disturbance in the dynamics.] we use the same terminal controllers as <cit.>, i.e., for i = 1,2,3:

    K_i^(1)   = K_i^(3) = [ -0.0544 -0.1398 ]
    
    		K_i^(2)   = K_i^(4) = [ -0.1544 -0.0295 ],

and the same robustly positive invariant terminal sets, i.e., for i = 1,2,3:

    𝒳_T, i = {x |    [  P; -P ]x ≤γ1, 
       P = [  7.8514  8.1971;  8.1957 -7.8503 ], γ=47}.

Finally, we compute local terminal costs V_f, i(x_i) = x_i^⊤Φ_i x_i, satisfying Assumption <ref>, using the procedure in Appendix <ref>, which gives

    Φ_1    = [ 12.67  8.87;  8.87  8.14 ]· 10^-4 + Φ
    Φ_2    = [ 10.58  7.90;  7.90  8.26 ]· 10^-4 + Φ
    Φ_3    = [ 8.53 5.43; 5.43 7.10 ]· 10^-4 + Φ
    Φ   = [ 3.938 1.262; 1.262 4.346 ].

For the SwA approach the horizon is N = 5 and the quadratic programs are solved using CasADi and the qrqp solver <cit.>.
Algorithm <ref> is run with T_ADMM=50 and ρ = 0.5.
These values are hand chosen to ensure the total residual is less than 0.01. 
For this system the edge case of continuous switching is never observed; T_cut is consequently set to also be 50.
For the initial guess in the first time step the trivial guesses u_i = (0,…,0) are used.
For the RCS approach the MIQPs are solved using Gurobi <cit.>.

Fig.  <ref> shows the trajectories under each approach for the initial condition, x_1 = [ -11 -18 ]^⊤, x_2 = [   2 -19 ]^⊤, x_3 = [ 15 19 ]^⊤.
Both approaches drive the states to the origin.
Fig.  <ref> shows, for 100 randomly sampled initial conditions, the performance difference between the two approaches

    Δ J = J_RCS - J_SwA

where J_RCS and J_SwA are the incurred closed-loop costs

    ∑_t = 0^30∑_i = 1^3 ℓ_i(x_i, u_i, {x_j}_j ∈𝒩_i).

The SwA approach generally performs slightly better, in terms of (<ref>), than the RCS approach, with occasional outlier initial conditions for which the performance improvement is much more significant and one initial condition in which the RCS approach performs better.
Fig.  <ref> also compares the online computation time required, averaged over all time steps in each simulation.
The SwA approach is much faster due to fact that only quadratic programs need to be solved, while the RCS approach involves solving MIQPs.
For all initial conditions the maximum computation time required by the SwA approach is less than the minimum computation time required by the RCS approach.


Now consider the same network, with the same terminal controllers, but with significantly higher coupling:

    A_12 = A_21 = A_23 = A_31 = [ 0.16    0;    0 0.16 ].

The uncontrolled global system is now unstable.
Under this coupling there do not exist local terminal sets that are robust to the coupling effects over the entire neighboring state spaces, and the RCS approach is hence not applicable.
There do, however, exist terminal sets satisfying Assumption <ref>.
Using the procedure in Appendix <ref>, initialized with the sets in (<ref>), we determine that the terminal sets in (<ref>) satisfy Assumption <ref> under the stronger coupling (<ref>).
Using the procedure in Appendix <ref> the terminal costs V_f, i(x_i) = x_i^⊤Φ_i x_i are computed as

    Φ_1    = [ 40.98 28.29; 28.29 43.73 ]Φ_2 = [ 32.07 20.90; 20.90 35.91 ]
             Φ_3 = [ 31.97 20.83; 20.83 35.07 ].

For the system with the increased coupling (<ref>) the continuous switching case is observed to occur rarely.
Algorithm <ref> is hence run with T_ADMM=75, T_cut=50, and ρ = 5, with these values again hand chosen to ensure that the total residual is less than 0.01.
Fig.  <ref> shows a closed-loop trajectory under this stronger coupling, demonstrating how our approach can provide stabilizing control for stronger coupling, due to its explicit handling of shared states and consequently less strict requirements on the terminal sets.
Furthermore, Fig. <ref> demonstrates the convergence of Algorithm <ref> in the third time step when the continuous switching edge case occurs.
Subsystem 1 continuously switches its local switching sequence, causing the local input and global residual to not converge.
For iterations larger than T_cut=50 the local switching sequence does not change, and the local input and global residual converge.






 §.§ Hybrid Vehicle Platooning

We now consider the problem of hybrid vehicle platooning that was presented as a benchmark problem for hybrid distributed MPC methods in <cit.>.
This example is not a simple regulation problem, and therefore the approach in <cit.> cannot be used.
Consider M vehicles that must form a platoon configuration while traveling in a single lane.
Each vehicle i has state x_i = [ p_i v_i ]^⊤, where p_i ∈ℝ and v_i ∈ℝ are the position and velocity of the vehicle, and control u_i ∈ℝ, the normalized throttle position.
The PWA vehicle dynamics x_i^+ = f_i(x_i, u_i) have seven PWA regions, modeling discrete gear changes that switch dependent on velocity, changing the traction force generated by the vehicle's control input.
For the sake of space the reader is referred to <cit.> for the full dynamics.
The control goal is for each vehicle to maintain a desired separation from the vehicle in front of it while the front vehicle follows a reference trajectory.
Without loss of generality it is assumed that the set of vehicles ℳ is ordered by the positions of vehicles, such that vehicle 1 is the front vehicle, and vehicle i tracks the trajectory of vehicle i-1.
The stage costs are then

    ℓ_1(x_1, u_1)    = x_i - r - η_Q + u_i_R 
    ℓ_i(x_i, u_i, x_i-1)    = x_i - x_i-1 - η_Q + u_i_R,

for i = 2,⋯,M with Q = [   1   0;   0 0.1 ], R = 1, η = [ 50  0 ]^⊤ and r ∈ℝ^2 the reference state.
Velocity and position bounds 

    3.94 ≤ v_i ≤ 45.84 
    
    		0 ≤ p_i ≤ 10000

define local convex constraints on states and inputs, and the coupling constraints

    h_i(x_i, x_i-1) = x_i - x_i-1 + d_safe≤ 0

enforce a safe distance d_safe = 25 between vehicles.
Vehicles are initialized with uniformly randomized velocities in the range [  5 30 ] and with uniformly randomized positions in the range [  50 100 ] behind the preceding vehicle. 

We compare our approach against existing controllers outlined in <cit.>.
For our approach, referred to as SwA, we again use a horizon of N = 5, T_ADMM = 100 and ρ = 0.5 for Algorithm <ref>, again ensuring the total residual is less than 0.01, and solve the quadratic programs with CasADi and the qpOASES solver <cit.>.
Once again the edge case of continuous switching is never observed; T_cut is consequently set to also be 100.
We implement a multi-start approach with two initial guesses, choosing the solution that yields a lower cost.
At time step t the shifted solution from time step t-1 is used, appending a constant control value,

    u_i = (u̅^⊤_i(1), …, u̅^⊤_i(N-1), u̅^⊤_i(N-1) )^⊤,

as well as a constant control sequence, maintaining the current velocity v_i.
The MPC controllers for comparison from <cit.> are based on converting the PWA dynamics to MLD form <cit.> and solving MIQPs.
The centralized controller, referred to as C, solves 𝒫_d directly as one MIQP, considering all vehicles in the one optimization problem, solving for all control inputs while considering all couplings.
The sequential controller, referred to as S, is a distributed controller where each vehicle solves a local MIQP once, in a predefined order, communicating the solutions down the sequence of vehicles.
Finally, the non-convex ADMM controller, referred to as NCA, applies ADMM directly to 𝒫_d, with vehicles iteratively solving local MIQPs.
For NCA 100 iterations are used for a fair comparison between the two ADMM-based approaches.
All MIQPs are solved with Gurobi <cit.>.

To further illustrate the mechanism of Algorithm <ref> Fig.  <ref> shows the local solutions u_i(0) throughout the first 20 iterations[For the remaining 80 iterations the values change negligibly.] at the first time step t = 0 of a simulation with M = 5.
Red dots indicate iterations in which local systems changed to new switching sequences, changing the global approximation 𝒫_s.
It is seen that the inputs vary sharply as vehicles iteratively solve their local minimization problems, communicate solutions, and identify new local switching sequences.
Eventually the changing of switching sequences stops, indicating that the local solutions are no longer pushed to the boundaries of the local constraint sets, and the solutions have settled in a local minimum of 𝒫_d.
Fig.  <ref> shows the total residual over the entire 100 iterations, demonstrating that subsystems come to agree on the values of the coupled states, as in Theorem <ref>.


Fig.  <ref> shows, for M = 10 vehicles, the cumulative tracking cost

    J(t) = ∑_τ = 0^t ℓ_1(x_1(t), u_1(t)) + ∑_i ∈ℳ∖{1}ℓ_i(x_i(t), u_i(t), x_i-1(t))
 
under each of the different controllers.
All controllers manage to reach a stable platoon formation, where the tracking cost no longer grows.
The centralized controller naturally performs the best as it finds the global optimum of 𝒫_d. 
The SwA controller outperforms both other distributed control solutions.
Fig.  <ref> shows the total tracking cost and the computation times for each controller with M = 5,10, and 15 vehicles.
As the number of vehicles increases the performance of the sequential and non-convex ADMM controllers decreases significantly with respect to the centralized controller, while the performance drop of the SwA controller scales much better.
The computation times show that the SwA is the fastest controller as it involves the solutions to only QPs rather than MIQPs.
Both the NCA and SwA approaches have computation times that scale elegantly with M, as the size of the optimization problems and the number of iterations is independent of M.
However, as NCA solves MIQPs iteratively rather than QPs, it has a significantly higher computational burden.





§ CONCLUSIONS

This paper has presented a new approach for distributed MPC for PWA systems.
The proposed controller is based on a novel switching ADMM procedure that solves a globally formulated non-convex optimal control problem, without requiring mixed-integer programming, and provably giving agreement on the values of shared states.
The recursive feasibility and stability of the resulting control scheme is proven under additional assumptions on the system.
The theoretical properties of the approach have been demonstrated on a small numerical example, whilst the practical utility has been demonstrated on a larger, more realistic, hybrid control problem.

Future work will look at quantifying the level of suboptimality between the closed loop performance of the proposed controller and the theoretically optimal performance of a centralized mixed-integer-based MPC controller.




§ COMPUTATION OF TERMINAL SETS IN ASSUMPTION <REF>


The procedure presented in <cit.> computes local positively invariant sets that are robust to the effects of coupling from neighboring sub-systems.
The procedure involves, for each sub-system i, computing a sequence of sets, beginning from a reasonably large polytope 𝒳_0, i⊆⋃_l ∈ℒ̂_i{x_i ∈ P_i^(l) | K_i^(l) x_i ∈𝒰_i}, terminating in the robustly positive invariant set 𝒳_T, i, which each set in the sequence being a contraction of the previous one.
See <cit.> for construction of 𝒳_0, i.

In this appendix a modified procedure with respect to that in <cit.> is presented for which the resulting robustly positive invariant sets are robust only to coupling for neighboring sub-systems within their own local terminal sets, as by Assumption <ref>.
A formal proof of the positive invariance, missing in <cit.>, is also provided.

Define the matrix A_cl, i^(l) = A_i^(l) + B_i^(l) K_i^(l). 
Define the set operation

    𝒬_i^(l)(𝒳, 𝒲)    = {x_i ∈ℝ^n_i | A_cl, i^(l) x_i + w ∈𝒳, ∀ w ∈𝒲}
       = (A_cl, i^(l))^-1 (𝒳⊖𝒲),

i.e., a one-step controllable set, robust to disturbances in the set 𝒲.
Algorithm <ref> outlines the procedure for computing 𝒳_T, i for i ∈ℳ.

[<cit.>]
	
	For all i ∈ℳ there exists τ∈ℕ^+ and τ < ∞ such that 𝒳_τ, i = 𝒳_τ-1, i and 𝒳_τ, i≠∅.


	Under Assumption <ref> the output sets from Algorithm <ref> are robustly positive invariant as in Assumption <ref>.


	As at each iteration τ of Algorithm <ref> the set 𝒳_τ, i shrinks, then for all i ∈ℳ, if x_i ∈𝒳_T, i and x_j ∈𝒳_T, j for all j ∈𝒩_i, then x_i ∈𝒳_τ, i and x_j ∈𝒳_τ-1, i for all j ∈𝒩_i and for all τ≥ 1.
	Then, by the definitions of 𝒬_i^(l) and 𝒲_τ-1, i, A_cl, i^(l) x_i + ∑_j ∈𝒩_i A_ij x_j ∈𝒳_τ-1, i for all l ∈ℒ_i and all τ≥ 2.
	Then A_cl, i^(l) x_i + ∑_j ∈𝒩_i A_ij x_j ∈𝒳_T, i for all l ∈ℒ_i.
	Hence, 𝒳_T, i is robustly positive invariant, robust to the coupling of neighboring sub-systems within their respective terminal sets, under arbitrary switching between the PWA regions l ∈ℒ_i.
	It follows directly that 𝒳_T, i is robustly positive invariant for the PWA sub-system under the true switching.



§ COMPUTATION OF V_F, I
 
Given the structures of ℓ_i and V_T, i (<ref>) can be written as

    ∑_i ∈ℳ   (⋆)^⊤Φ_i (  (A_i^(l_i) + B_i^(l_i) K_i^(l_i)) x_i + ∑_j ∈𝒩_i A_ij^(l_i)x_j ) 
       - x_i^⊤Φ_i x_i + x_i^⊤ Q_i x_i + (⋆)^⊤ R_i (K_i^(l_i)x_i) 
       + ∑_j ∈𝒩_i x_j^⊤ Q_ij x_j ≤ 0, x_i ∈ P_i^(l_i),

for all x_i ∈𝒳_T, i, l_i ∈ℒ̂_i, i ∈ℳ, recalling that c_i^(l_i) = 0 for l_i ∈ℒ̂_i.
Finding Φ_i's such that (<ref>) is satisfied can be achieved by considering the system from a global view, ensuring the inequality holds for every combination of PWA regions of the subsystems.
This is equivalent to considering the global system as a single large PWA system with ∏_i ∈ℳL_i regions and using the approach presented in <cit.> for single PWA systems.
Condition (<ref>) can be reformulated as the set of LMIs

    A_cl^⊤(l)    Φ A_cl(l) - Φ + Q + K(l) R K(l) ≺ 0, 
       ∀ l = (l_1,…,l_M) ∈ℒ̂_1 ×…×ℒ̂_M

where Φ = blck(Φ_1, …, Φ_M), Q = blck(Q_1 + ∑_j|1 ∈𝒩_j Q_j1,…,Q_M + ∑_j | M ∈𝒩_jQ_jM), R = blck(R_1, …, R_M), K(l) = blck(K_1^(l_1), …, K_M^(M)), and

    A_cl(l)    = [  A_cl, 1^(l_1)     A_12^(l_1)              …     A_1M^(l_1);     A_21^(l_2)              ⋱                             ⋮;              ⋮                             ⋱               ;     A_M1^(l_M)                               A_cl, M^(l_M). ]

with A_ij^(l_i) = 0 if j ∉𝒩_i.
Evidently, (<ref>) requires centralized computation, and may become intractable for large M or |ℒ̂_i|.
We leave, however, distributed and consequently less computationally demanding solutions to future work.
Finally, (<ref>) is linear in the matrix Φ as the linear controllers K_i^(l) have been assumed to be computed a priori as in <cit.>.
The linear controllers could be computed simultaneously using the same expression, leveraging the Schur complement and a change of variables to maintain an LMI expression <cit.>.



§ ACKNOWLEDGMENT

The authors would like to thank the authors of <cit.>, A. Ma, D. Li, and Y. Li, for helpfully providing the code for their approach.


IEEEtran




[
    < g r a p h i c s >
]Samuel Mallick 
	received the B.Sc. and M.Sc. degrees from The University of Melbourne in 2020 and 2022, respectively. He is currently a Ph.D. candidate at the Delft Center for Systems and Control, Delft University of Technology, The Netherlands.
	
	His research interests include model predictive control, reinforcement learning, and distributed control of large-scale and hybrid systems.





[
    < g r a p h i c s >
]Azita Dabiri 
	received the Ph.D. degree from the Automatic Control Group, Chalmers University of Technology, in 2016. She was a Post-Doctoral Researcher with the Department of Transport and Planning, Delft University of Technology, from 2017 to 2019. In 2019, she received an ERCIM Fellowship and also a Marie Curie Individual Fellowship, which allowed her to perform research at the Norwegian University of Technology (NTNU), as a Post-Doctoral Researcher, from 2019 to 2020, before joining the Delft Center for Systems and Control, Delft University of Technology, as an Assistant Professor. Her research interests are in the areas of integration of model-based and learning-based control and its applications in transportation networks.




[
    < g r a p h i c s >
]Bart De Schutter(Fellow, IEEE) 
	received the PhD degree (summa cum laude) in applied sciences from KU Leuven, Belgium, in 1996. He is currently a Full Professor and Head of Department at the Delft Center for Systems and Control, Delft
	University of Technology, The Netherlands. His research interests include multi-level
	and multi-agent control, model predictive control, learning-based control, and control
	of hybrid systems, with applications in intelligent transportation systems and smart energy systems. 
	
	Prof. De Schutter is a Senior Editor of the IEEE Transactions on Intelligent Transportation Systems.




