ADMM for MPC with state and input constraints, and input nonlinearity
Sebastian East Mark Cannon

 On Behind the Physics of the Thermoelectricity of Topological Insulators 
    Daniel Baldomir^1∗ and Daniel Faílde^1



^1Departamento de Física Aplicada, Instituto de Investigacións Tecnolóxicas. Universidade de 

Santiago de Compostela E-15782 Campus Vida s/n, Santiago de Compostela, Spain



^∗daniel.baldomir@usc.es

    
=======================================================================================================================================================================================================================================================

empty
empty




In this paper we propose an Alternating Direction Method of Multipliers (ADMM) algorithm for solving a Model Predictive Control (MPC) optimization problem, in which the system has state and input constraints and a nonlinear input map. The resulting optimization is nonconvex, and we provide a proof of convergence to a point satisfying necessary conditions for optimality. This general method is proposed as a solution for blended mode control of hybrid electric vehicles, to allow optimization in real time. To demonstrate the properties of the algorithm we conduct numerical experiments on randomly generated problems, and show that the algorithm is effective for achieving an approximate solution, but has limitations when an exact solution is required.







§ INTRODUCTION


In Model Predictive Control (MPC), control inputs are optimized by solving an open-loop optimal control problem over a finite prediction horizon, whilst explicitly accounting for constraints on system states and inputs. This approach is widely used in process control applications, in which slow plant dynamics allow the optimization to be solved in real time, and the control of mechanical and electrical systems is becoming more feasible with improvements in embedded controller hardware <cit.>. To further increase MPC applicability, recent research has considered exploiting properties of optimization problems to improve computational efficiency.

Active set <cit.> and interior point methods <cit.> are commonly employed in MPC, but for reduced computational complexity first order approaches such as fast gradient methods <cit.> and alternating direction method of multipliers (ADMM) <cit.> have received recent attention. ADMM has been shown to be particularly effective in distributed applications (D-ADMM) <cit.> where the separability of the cost function can be leveraged. This paper likewise exploits separability of the cost and constraints, although distributed optimization and control are not considered here. We refer the reader to <cit.> for a full description of ADMM and a survey of applications.

This paper is motivated by the problem of minimizing fuel consumption in a hybrid electric vehicle, achieved by distributing the demanded load between two available power sources (an internal combustion engine and an electric motor) in a blended mode <cit.>. The problem is constrained by limits on available power and bounds on the battery state of charge, and MPC is a suitable framework since its feedback mechanism provides a degree of robustness to discrepancies between the predicted and actual power demanded. Although dynamic programming is typically used for the solution of the corresponding optimization <cit.>, the computation required for a sufficiently accurate solution prohibits a real-time implementation. A projected Newton method for solving this problem subject to a terminal state of charge constraint was proposed in <cit.>, but the method becomes intractable for the general case of constraints on the state charge at each instant on a prediction horizon. Here, we use ADMM to leverage a separable cost function to enforce constraints at all  future time-steps.

The paper is organised as follows: Section 2 gives the problem definition; Section 3 states the proposed ADMM algorithm; Section 4 describes the energy management problem more fully and discusses its solution for randomly generated examples; Section 5 provides conclusions. Convergence and optimality analyses are included in the appendix.
 


§ PROBLEM DEFINITION


Consider a system represented by the discrete time model

    x_k+1 = A_kx_k+b_k(u_k)

where x_k ∈ℝ^n_x is the system state, u_k ∈ℝ^n_u is the control input, A_k∈ℝ^n_x× n_x, and b_k(·) is convex. The states and control inputs are subject to elementwise bounds and the associated optimal control problem is defined 

    x,uminimize       f(x) + g(u)
       subject to       x_k+1 = A_kx_k+b_k(u_k)  
    x_k+1≤ x_k+1≤x_k+1
    u_k ≤ u_k ≤u_k
      
    k = 0,…,N-1

where x=(x_1,…,x_N) is the vector of future states over a prediction horizon of N steps, u=(u_0,…,u_N-1) is the vector of control inputs over the prediction horizon, and f,g are state and control cost functions. It is assumed that f and g are convex and separable. We say that a function a(c), where c = (c_1,…,c_N_c), is separable if

    a(c) = ∑_i=1^N_c a_i(c_i) .


Note that problem (<ref>) is in general nonconvex due to the presence of nonlinear dynamics and state constraints.


Equality constraints: x and x_k can be written in terms of
b(u)=[b_0^⊤(u_0) ⋯ b_N-1^⊤(u_N-1)]^⊤ as

    x    = Φ x_0 + Ψ b(u) 
    
    x_k   = Φ_k x_0 + Ψ_k b(u)

 where Φ_k and Ψ_k are the kth block rows of Φ and Ψ:

    Φ = [          A_0;            ⋮; ∏_k=0^N-1A_k ] ,
      Ψ = [            I            …            0;            ⋮            ⋱            ⋮; ∏_k=1^N-1A_k            …            I ] .

Here Ψ∈ℝ^N n_x × N n_x is lower-triangular with ijth block given by ∏_k=j^i-1A_k for j<i and by the identity matrix for i=j, and the ith block-row of Φ∈ℝ^N n_x × n_x is ∏_k=0^i-1A_k.


























Inequality constraints:
Define elementwise indicator functions h_k^x,  h_k^u and h^x, h^u as

h_k^x(x_k)    =  0     x_k ∈[x_k,x_k] 
 ∞    otherwise 
    h^x(x)   =∑_k=1^Nh_k^x(x_k) , 

h_k^u(u_k)    =  0     u_k ∈[u_k,u_k] 
 ∞    otherwise 
    h^u(u)    = ∑_k=0^N-1h_k^u(u_k) .

Using these definitions we rewrite (<ref>) in a more convenient form without explicit inequality constraints as 

    x,uminimize       f(x) + g(u) + h^x(x) + h^u(u) 
       subject to       x = Φ x_0 + Ψ b(u) .






























§ ADMM ALGORITHM

Introducing v as a substitute for b(u), we rewrite (<ref>) as

    u,x,vminimize       f(x) + g(u) + h^x(x) + h^u(u)
       subject to      Φ x_0 + Ψ v - x = 0 
              b(u) - v = 0 .

The associated augmented Lagrangian function is 

    L(u,v,x,y,z) = f(x) + g(u) + h^x(x) + h^u(u) 
    
    +ρ_1/2b(u) - v +y ^2 + ρ_2/2Φ x_0 + Ψ v - x +z ^2 ,

and the ADMM iteration is obtained (see e.g. <cit.>) as

   u_k^j+1        =  u_kmin( h_k^u(u_k) + g_k(u_k) + ρ_1/2 (b_k(u_k) - v^j_k + y^j_k )^2 ) 

   v^j+1        = [t] v min( ρ_1/2 b(u^j+1)    - v +y^j ^2 
 
    + ρ_2/2Φx_0 + Ψv - x^j +z^j ^2 )  

   x_k+1^j+1        = [t] x_k+1 min(
    h_k^x(x_k+1) + f_k+1(x_k+1) 

     +ρ_2/2(Φ_k+1 x_0 + Ψ_k+1 v^j+1- x_k+1 +  z^j_k )^2 )  

   y^j+1        = y^j + b(u^j+1) - v^j+1 

   z^j+1        = z^j + Φx_0 + Ψv^j+1-x^j+1

for k = 0,…,N-1, where j is the iteration counter. We assume that a suitable solver is available for updating u^j+1 and x^j+1. For the application considered in Section <ref> analytical solutions exist for the minimizers u_k^j+1 and x_k+1^j+1, which are computed by finding the roots of cubic equations and projecting these onto the inequality constraints in (<ref>).



The update for v can be written explicitly as

    v^j+1 =  ( ρ_1I+ρ_2 Ψ^⊤Ψ)^-1[    ρ_1(b(u^j+1) + y^j) 
       + ρ_2Ψ^⊤(-Φ x_0 +x^j-z^j)] .

For the application considered in Section <ref>, where A_k = 1, the matrix Ψ is a lower triangular matrix of 1s, so (ρ_1I+ρ_2Ψ^⊤Ψ)^-1 becomes a linear time-invariant filter. Approximating this matrix by setting the elements that do not exceed a given threshold to zero then yields a banded matrix, making the computation of v^j+1 straightforward. For the general case in which the matrices A_k defining Ψ are arbitrary, the results in <cit.> and <cit.> provide conditions under which (ρ_1 I + ρ_2 Ψ^⊤Ψ)^-1 is approximately banded and bound the rates of decay of elements with distance from the diagonal. 


The iteration is terminated when primal and dual residual variables have fallen below pre-defined thresholds,
ϵ^primal and ϵ^dual,
chosen based on the required accuracy of the optimization and the typical magnitudes of decision variables:

    r^j+1_2 ≤ϵ^primal ,   s^j+1_2 ≤ϵ^dual .







In the appendix we analyse the optimality and convergence properties of the algorithm, and provide definitions of the residuals r^j and s^j. 

We also give conditions for the algorithm to converge to the global minimum, despite the problem being nonconvex, however this condition cannot generally be determined a priori, and the algorithm may converge to a local minimum, or even a maximum, if it is not met.




§ APPLICATION TO PHEV SUPERVISORY CONTROL

We are interested in the problem of minimizing the fuel consumption of a plug-in hybrid electric vehicle (PHEV) by controlling the power balance between the electric motor and internal combustion engine in a blended mode. The power, d, demanded by the driver is split between the motor m and engine u so that u+m = d. The system can be modelled as 

    x_k+1 = x_k+b_k(u_k)

where x_k is the state of charge of the battery, b_k accounts for the charge  loss-dynamics, and the engine power u_k is the control input (i.e. n_x=n_u = 1 and A_k = 1). The fuel minimization problem over a horizon of N time steps is 

    uminimize      ∑_k=0^N-1 g_k(u_k) 
       subject to       x_k+1 = x_k+b_k(u_k) 
    x_k+1≤ x_k+1≤x_k+1
    u_k ≤ u_k ≤u_k
     k = 0,…,N-1

where g_k is the fuel consumption determined from a quasi-static model at time-step k; u_k and u_k represent upper and lower limits on the engine power  output; and x_k+1 and x_k+1 represent the bounds on state of charge. The cost and input loss functions are represented by time-varying quadratic functions as

    g_k(u_k)    = α_2,ku_k^2 + α_1,ku_k + α_0,k
    
    b_k(u_k)    = -β_2,k(d_k - u_k)^2 - β_1,k(d_k - u_k) - β_0,k

See <cit.> for details on how these functions are obtained from the fuel map and electrical loss maps for a particular PHEV.

To demonstrate the performance of the ADMM algorithm without reference to a particular PHEV powertrain, we generate random sets of nominal systems of the same form as the motivating problem in (<ref>). We generate uniformly distributed random disturbance values d_k ∈ [-1,1]; uniformly distributed α and β coefficients α_2,k,β_2,k∈ [0,0.1], α_1,k,β_1,k∈ [0,1], and α_0,k,β_0,k =0, and a uniformly distributed initial state x_0 ∈ [-0.5,0.5]. We use u_k =-0.5 and u_k =0.5 as control input limits, and state limits of x_k = -2 and x_k = 2. The ADMM iteration is initialised assuming that the inequality constraints in (<ref>) are inactive:



    u_k^0    =  π_k^u[ u_kmin  g_k(u_k) ], 
          v^0     = b(u^0),
    
    
    x_k+1^0    = π _k+1 ^x [  Φ_k+1 x_0 + Ψ_k+1 v^0   ],
          y^0     = 0,
       
    z^0 = 0,










for k = 0,…,N-1, 
where π_k^u, π_k^x are projection operators:

    π_k^u(u_k)    = min{u_k,max{u_k,u_k }}
    
    π_k^x(x_k)    =min{x_k ,  max{x_k, x_k}} .


If Φ x_0 + Ψ v^0 = x^0, then these initial values are optimal and the state constraints are inactive. The stopping conditions used are equal i.e ϵ^primal=ϵ^dual=ϵ. The simulations are run in Matlab, using an Intel 2.60 GHz i7-6700HQ CPU.



 §.§ Variation of iteration count with ρ

The rate of convergence of the ADMM algorithm is affected by the values of ρ_1 and ρ_2. For the algorithm to be practically useful, it is desirable for the optimal values for these parameters to be independent of both horizon length and system parameters, so that a single set of values can be used in all conditions. Here we demonstrate the effect of ρ_1 and ρ_2 on the convergence properties of the algorithm for 20 randomly generated systems. Figure <ref> shows the average numbers of iterations required for completion of the algorithm for horizon lengths of 50, 100, 200, and 400 with ϵ = 0.001, 10^-1≤ρ_1≤ 10^2 and 10^-2≤ρ_2≤ 10.




For each horizon length, the minimum number of iterations for convergence is obtained with 1 ≤ρ_1 ≤ 10 and 0.1 ≤ρ_2 ≤ 1. Moreover, variations in ρ within an order of magnitude of this minimum do not produce large increases in total number of iterations, implying a degree of robustness. Therefore ρ values within these ranges provide close to the minimum number of iterations for horizon lengths between 50 and 400, suggesting that a single set of values is appropriate for the PHEV problem.

Although the optimal ρ values are similar, the minimum number of iterations for completion increases significantly with horizon length. This is explored further in the following sections.





 §.§ Variation of required number of iterations with ϵ

The threshold ϵ determines the degree of optimality at termination, since (u,x) → (u^*,x^*) as ϵ→ 0 (see Appendix A). However, a trade-off exists between accuracy and computation since small values of ϵ require many iterations. Here we investigate how the number of iterations before termination varies with ϵ for 200 randomly generated problems with N=100. Figure <ref> shows the number of iterations required with 10^-4≤ϵ≤ 10^0, ρ_1=1, and ρ_2=0.2.

Clearly, for a low-accuracy solution only a few tens of iterations are required, but this rises rapidly for high-accuracy solutions. For example, with ϵ =0.14  the median iteration count is 29, but ϵ = 10^-4 requires a median iteration count of 210. Furthermore, the increase in iteration count is highly nonlinear as ϵ is reduced, and for small ϵ modest further reductions have a large effect on the iteration count. 

The variation in the number of iterations required increases as ϵ is reduced (Figure <ref>): the 98th percentile is 45 for ϵ =0.14 (∼ 2 times the median), rising to 2374 (∼ 10 times the median) for ϵ = 10^-4. Again, this increase in uncertainty increases nonlinearly as ϵ is reduced. Furthermore, the distribution of the data is heavily skewed, and outliers become more extreme as ϵ is reduced. This is problematic as the MPC optimization algorithm must be designed for the worst case, and at low values of ϵ such a high potential iteration count could be extremely prohibitive.

These properties suggest that the ADMM algorithm may be appropriate for finding an approximate solution to the PHEV problem that is then used to initialise another method, but inappropriate for obtaining an exact solution. For example, the ADMM algorithm could be used to obtain the active set, then an alternative method (such as that proposed in <cit.>) used to obtain the corresponding optimum. This observation is application-specific as large numbers of iterations may not be problematic if the system dynamics are sufficiently slow or computational resources are sufficiently high.







 §.§ Variation of computation with horizon length

The MPC framework can be implemented in the PHEV energy management problem with either a receding horizon (with fixed N) or a shrinking horizon (with N reducing as the vehicle progresses through a journey). Although the algorithm can be tuned for a specific horizon length in the receding horizon case, for a shrinking horizon implementation it must be robust to a range of horizon lengths. Thus, for a 1 hour journey sampled at 1 s, N varies between 1 and 3600. Here we demonstrate the effect of horizon length on number of iterations and time to termination for 200 systems, with N varying between 10 and 400. We choose ϵ=10^-2, and ρ_1=1, ρ_2=0.2. The results are shown in Figure <ref>.

For small N, increases in horizon length produce large increases in the number of iterations (the median for N=25 is 22, increasing to 40 when N = 75), but this effect tapers off quickly, and the median iteration count is 105 for N = 350. There is no clear dependence on N of the uncertainty in the required number of iterations; the 98th percentile has a maximum value of 273 iterations at N=80, a minimum of 85 at N=10, and large variations between these values.

Although the iteration count does not increase significantly with longer horizons, the total time taken to find the solution increases roughly linearly with horizon length, from a median of 0.0075 s for N=10, to a median of 1.39 s for N=400. The uncertainty also appears to increase within a near-linear envelope. This is due to the computation required to update u^j+1, for which N cubic equations must be solved at each iteration. This suggests that for long prediction horizons, improving the solution times for u^j+1 and x^j+1 would have a significant effect on total solution time.


Note that the absolute computation times presented here demonstrate the relationship between horizon length and computation, but 

the hardware used for the simulations does not correlate to that typically found in a PHEV, and Matlab is not typically used for embedded computing.







§ CONCLUSIONS



A general MPC optimization problem is considered for systems with separable cost functions, nonlinear input map and simple constraints, motivated by the blended mode energy management problem for PHEVs. We propose an ADMM algorithm for solving this problem and demonstrate conditions for convergence to a solution satisfying first order necessary conditions for optimality.

Numerical experiments show that an approximate solution is reached within a few tens of iterations, while significantly more are required for an exact solution. The effect of algorithm parameters on the number of iterations and solution times is also investigated.

Future work will investigate methods for implementing the proposed algorithm in a PHEV, to obtain globally optimal control inputs that are implementable in real time.  





§ APPENDIX




 §.§ Optimality

We use a similar method to that given in Section 3 of <cit.>, however we make use of an alternative Lyapunov function to demonstrate convergence, as required by the nonconvexity of the problem. Problem (<ref>) is first rewritten as 

    uminimize      f̂(x̂) + ĝ(u)  
       subject to      b̂(u) + Bx̂ = 0

where

    f̂(x̂) = f(x) + h^x(x),  ĝ(u) = g(u) + h^u(u)
       b̂(u) = [ [  b(u); Φ x_0 ]],  
    B = [ [ -I  0;  Ψ -I ]], x̂ = [ [ v; x ]] .

Then, by expressing the augmented Lagrangian (<ref>) as

    L(u,x̂,ŷ) = f̂(x̂) + ĝ(u) + ŷ^⊤(b̂(u)+Bx̂)  + 12b̂(u)+Bx̂^2_R

where w ^2_R = w^⊤ R w for any vector w of conformal dimensions, and defining ŷ = [y^⊤  z^⊤]^⊤, R = {ρ_1I,ρ_2 I},





the ADMM iteration becomes


    u^j+1   = umin{ĝ(u) + (ŷ^j)^⊤(b̂(u) + Bx̂^j)  
        + 12b̂(u)+Bx̂^j ^2_R }
    x̂^j+1   = x̂min{f̂(x̂) + (ŷ^j)^⊤(b̂(u^j+1) + Bx̂)  
        + 12b̂(u^j+1)+Bx̂^2_R }
    ŷ^j+1    =  
                      ŷ^j + R(b̂(u^j+1)+Bx̂^j+1)



Slater's condition holds, i.e. u,v,x exist such that
u_k  < u_k < u_k, x_k+1 < x_k+1 < x_k+1 for all k, where x = Φ x_0 + Ψ v and v = b(u).



The first order conditions satisfied by a (possibly local, possibly maximal) optimal solution (u^*,x̂^*,ŷ^*) are
 

    0 =b̂(u^*) + Bx̂^* 
        0 = ∇ĝ(u^*) + [∂b̂(u^*)]^⊤ŷ^* 
        0 = ∇f̂(x̂^*) + B^⊤ŷ^*


where ∇ĝ, ∇f̂ are the gradients of ĝ, f̂, and ∂b̂ is the Jacobian matrix of b̂.


(u^j+1,x̂^j+1,ŷ^j+1) satisfies the first order necessary  optimality conditions (<ref>) iff r^j+1=0 and s^j+1=0.



First consider condition (<ref>). Since x̂ = x̂^j+1 is by definition the minimizer of L(u^j+1,x̂,ŷ^j), we have

    ∇f̂(x̂^j+1) + B^⊤ + B^⊤ R (b̂(u^j+1) + B  ) = 0

and hence the update law =  + R(b̂(u^j+1) + B) implies ∇f̂() + B^⊤ = 0 for all j. Next, consider (<ref>). By definition u = u^j+1 minimizes L(u,,) and therefore u^j+1 necessarily satisfies

    ∇ĝ(u^j+1) + [∂b̂(u^j+1) ]^⊤(  + R (b̂(u^j+1) + B)) = 0

Hence the update law for  gives

    ∇ĝ(u^j+1) + [∂b̂(u^j+1) ]^⊤(  + RB ( - ) = 0

and we have, for all j,

    ∇ĝ(u^j+1) + [∂b̂ (u^j+1)]^⊤ + s^j+1 = 0

Therefore, the ADMM iteration converges to a point satisfying the first order necessary conditions (<ref>)-(<ref>) if and only if the residual variables defined by

    r^j+1   = b̂(u^j+1) + B 
    
    s^j+1   = [∂b̂(u^j+1) ]^⊤ RB ( - )

converge to zero.




 §.§ Convergence


Let (u^Ø,x̂^Ø,ŷ^Ø) denote a solution of (<ref>) that achieves the globally minimum value of the objective. Let L_0(u,,) = f̂() + ĝ(u) + ^⊤ (b̂(u) + B  ) and define _c as the set

    _c ={ (u,,) :  - ^Ø^2_R^-1 + B( - ^Ø) ^2_R + r(u,)^2_R ≤ c }

for c>0, where r(u,) = b̂(u) + B.

_c, (u^Ø,^Ø,^Ø), (u^j,^j,^j) are such that:


  * iteration (<ref>) is initialised at (u^0,^0,^0) ∈_c, i.e.
 
c ≥^0 - ^Ø^2_R^-1 + B(^0 - ^Ø^2_R + r(u^0,^0)^2_R,

  * u = u^j+1 in (<ref>a) is the global minimiser of L(u,^j,^j).



 
(u^j,^j,^j) ∈_c,
 for all j ∈ℤ^+, and r^j → 0 and s^j → 0 as j→∞.


We first show that V^j = V(u^j,^j,^j) is a Lyapunov function, where V is defined by

    V(u,,) =  - ^Ø^2_R^-1  + B( - ^Ø )^2_R + r(u,)^2_R .

This part of the proof is split into three steps:


  * The update law for u^j+1 can be equivalently written as

    u^j+1 = umin{ĝ(u) + ( + RB(-))^⊤b̂(u) }

and, from Assumption 3(b) it follows that 

    ĝ(u^j+1)    + (ŷ^j+1 + RB( - ))^⊤b̂(u^j+1) 
       ≤ĝ(u^Ø) + (ŷ^j+1 + RB( - ))^⊤b̂(u^Ø) .

Similarly, ∇f̂() + B^⊤ = 0 implies that

    f̂( ) + ( )^⊤ B ≤f̂(x̂^Ø ) + ( )^⊤ B x̂^Ø .

Therefore, defining p^j+1 = f̂ (  ) + ĝ(u^j+1) and p^Ø = f̂(x̂^Ø) + ĝ(u^Ø), we obtain

    p^j+1 - p^Ø = f̂() - f̂(x̂^Ø) + ĝ(u^j+1) - ĝ(u^Ø) 
       ≤ ()^⊤ (B(x̂^Ø - ) + b̂(u^Ø) - b̂(u^j+1)) 
          + ( -  )^⊤ B^⊤ R (b̂(u^Ø ) - b̂(u^j+1)) 
       = -()^⊤ r^j+1
          + ( - )^⊤ B^⊤ R (-r^j+1 + B( - x̂^Ø))



  * Since (u^Ø,x̂^Ø,ŷ^Ø) must satisfy the first order necessary conditions (<ref>), it follows that u=u^Ø is the minimiser of ĝ(u) + (ŷ^Ø)^⊤b̂(u) and = ^Ø is the minimiser of f̂( ) + ( ^Ø)^⊤ B. Therefore, we have, for all u and ,

    ĝ(u^Ø ) + (^Ø)^⊤b̂(u^Ø )    ≤ĝ(u) + (^Ø )^⊤b̂ (u) 
    f̂(^Ø ) + (^Ø )^⊤ B ^Ø   ≤f̂(x̂) + (^Ø )^⊤ B  ,
 
so that f̂(^Ø ) + ĝ ( u^Ø) ≤f̂ ( ) + ĝ (u) + (^Ø )^⊤ (B + b̂(u)) and hence p^Ø - p^j+1≤ (^Ø )^⊤ r^j+1.


  * Combining the bounds on p^Ø - p^j+1 in i) and ii) yields

    (^j+1 - ^Ø )^⊤ r^j+1 - ( -  )^⊤ B^⊤ R r^j+1
    
    + ( -  )^⊤ B^⊤ R B ( - ^Ø) ≤ 0 .

The first term in (<ref>) can be simplified using the update law =  + R r^j+1 and completing the square:

    2( - ^Ø )^⊤ r^j+1 = 2( -^Ø)^⊤ r^j+1 + 2  r^j+1^2_R 
       =  - ^Ø + R r^j+1 ^2 _R^-1 -  - ^Ø^2_R^-1 +  r^j+1^2_R 
       =  - ^Ø + ^j+1 - ^j  ^2 _R^-1 -  - ^Ø^2_R^-1 + r^j+1^2_R 
       = ^j+1 - ^Ø^2_R^-1 - ^j - ^Ø^2_R^-1 + r^j+1^2_R .

To simplify the 2nd and 3rd terms in (<ref>) we use

    r^j+1^2_R - 2( - )^⊤ B^⊤ R r^j+1
            + 2( -  )^⊤ B^⊤ R B ( - ^Ø) 
       = r^j+1 - B( - ) ^2_R  
            +  B( - ^Ø ) ^2_R - B( - ^Ø ) ^2_R .

Therefore (<ref>) is equivalent to

    ^j+1 - ^Ø^2_R^-1 - ^j - ^Ø^2_R^-1
    
    + r^j+1 - B(^j+1 - ^j) ^2_R + B( - ^Ø ) ^2_R 
    
     -  B( - ^Ø ) ^2_R ≤ 0

which, using the definition of V^j = V(u^j,^j,), is equivalent to 

    V^j+1 - V^j    ≤ - r^j+1 - B( - ) ^2_R 
          - r^j ^2_R + r^j+1^2_R 
        =  -r^j ^2_R - B( - ) ^2_R 
          +2(r^j+1)^⊤ R B ( - ) 
        = - r^j ^2_R -  B( - ) ^2_R 
          + 2( - )^⊤ B ( - ) .

To bound the RHS of this expression, we recall from i) that = and = are the minimizers of f̂() + ()^⊤ B and f̂() + ()^⊤ B, so that

    f̂() + ()^⊤ B ≤f̂ (  ) + ()^⊤ B

and

    f̂() + ()^⊤ B ≤f̂ (  ) + ()^⊤ B

By summing these two inequalities we obtain ( - )^⊤ B ( - ) ≤ 0, and it follows that 

    V^j+1 - V^j ≤ - r^j ^2_R -  B( -  ) ^2_R .


The positive invariance of _c follows from its definition as a level set of V(u,,) and from the fact that V^j is monotonically non-increasing with j. Asymptotic convergence of the residuals r^j and s^j can be established by summing both sides of the inequality satisfied by V^j+1 - V^j over all j ≥ 0:

    ∑_j=0^∞( r^j^2_R + B( - ) ^2_R ) ≤ V^0.

From this bound we can conclude that, as j →∞, r^j → 0 and B( - ) → 0, and hence also s^j → 0.


Lemma <ref> and Theorem <ref> imply that the ADMM iteration converges to a point satisfying the first order necessary optimality conditions, and this will be the global minimum if the problem is convex. The invariant nature of _c also implies that even when the problem is nonconvex, the algorithm will converge to the global minimum if it is initialised within a sub-level set of V that contains no local minima or maxima. This condition is, however, impossible to demonstrate a priori, and if the algorithm is initialised outside of this set it could converge to a local minimum or even a local maximum point. We note finally that the Theorem <ref> also demonstrates that the iteration necessarily terminates after a finite number of steps whenever the threshold ϵ on r^j+1 and s^j+1 is set to a fixed non-zero value.










unsrt



