Matrix Completion Using Alternating Minimization for Distribution System State Estimation
Yajing Liu April Sagan Andrey Bernstein Rui Yang Xinyang Zhou Yingchen Zhang

Brain Tumor Synthetic Segmentation in 3D Multimodal MRI Scans
    Mohammad Hamghalam1,20000-0003-2543-0712 Baiying Lei1 Tianfu Wang1
    July 2019
======================================================================


This paper examines the problem of state estimation in power distribution systems under low-observability conditions. The recently proposed constrained matrix completion method which combines the standard matrix completion method and power flow constraints has been shown to be effective in estimating voltage phasors under low-observability conditions using  single-snapshot information.  However, the method requires solving a semidefinite programming (SDP) problem, which becomes computationally infeasible for large systems and if multiple-snapshot (time-series) information is used. This paper proposes an efficient algorithm to solve the constrained matrix completion problem with time-series data. This algorithm is based on reformulating the matrix completion problem as a bilinear (non-convex) optimization problem, and applying the alternating minimization algorithm to solve this problem. This paper proves the summable convergence  of the proposed algorithm, and demonstrates its efficacy and scalability via IEEE 123-bus system and a real utility feeder system. This paper also explores the value of adding more data from the history in terms of computation time and estimation accuracy. 












§ INTRODUCTION

State estimation in power systems is critical for maintaining normal and secure operating conditions. In power transmission systems, the weighted least squares (WLS) method <cit.> is  well developed and widely used by utilities for state estimation; however, the limited availability of measurements renders the WLS inapplicable in  power distribution systems. The requirement for state estimation in  distribution systems is becoming stringent because the increased penetration of distributed energy resources—such as solar photovoltaics, wind turbines, and energy storage systems—on distribution systems introduces bidirectional power flows that might impact the system responses to various types of disturbances.

Much work has been done in an attempt to address the low observability state estimation problem. Reference <cit.> determined the optimal minimum measurement locations and used the difference between the measured and calculated voltages and complex powers to obtain the voltage profile of the whole network. Reference <cit.> used  the voltage sensitivity at each bus to the load at all buses to determine the voltage given the measured load at each location. References <cit.> applied neural networks, which do not require system models, to estimate the voltage in the distribution grid; however, these approaches still require the installation of significant numbers of phasor measurement units (PMUs) when considering large systems. References <cit.> used pseudo-measurements, which are estimates based on historical data to overcome the requirement of the installation of more meters; however, it is known that the pseudo-measurements typically have  larger measurement errors than the metered measurements.

An alternative approach that does not require installing new sensors or explicit computation of pseudo-measurements was recently proposed by <cit.>. This method leverages the standard matrix completion method <cit.> and augments it with linearized power flow constraints to acknowledge for the physical network constraints.
This method was shown to be effective in estimating voltage phasors using single-snapshot information.  However, this method requires solving a semidefinite programming (SDP) problem, which becomes computationally infeasible for large systems or if multiple-snapshot (time-series) information is used.





In fact,  temporal correlation (dependency between measurements at different time steps) exists among data in addition to the spatial correlation (dependency between measurements at different locations) and the correlation among measurement types (characterized by  power flow equations). By modeling both the temporal and spatial correlation of the data, <cit.> apply the standard matrix completion method to recover missing PMU data,  and <cit.> applies the standard matrix completion method with Bayesian estimation to recover missing low-voltage distribution system data.





This paper proposes an efficient algorithm to solve the constrained matrix completion problem with time-series data. To this end, this paper formulates  matrix X using time-series data  by considering both the  spatial  and  temporal correlation of the data. Similar to <cit.>, this paper leverages a linear approximation of the power flow equations as constraints to characterize the correlation among different measurement types.

Then, this paper applies 
the alternating minimization method,  which was a winner in the Netflix Challenge <cit.>,  to solve the proposed matrix completion problem for state estimation in multiphase distribution systems.

To apply the alternating minimization method,  matrix X is written in a bilinear form, i.e., X=UV, and then the algorithm finds the best U and V in an alternating fashion. Because of the bilinear term in the objective function, the problem becomes nonconvex with respect to any variable and thus the alternating minimization method is not guaranteed to converge to the global optima. However, 
results on the  convergence of the method were established under some regularity conditions in <cit.>. This paper proves that the objective function satisfies these conditions, which guarantee the convergence of the proposed method. 







The rest of this paper is organized as follows. Section <ref> introduces the standard matrix completion model, the formulation of the data matrix using time-series data, the linearized power flow constraints, and the full matrix completion model for state estimation in multiphase distribution systems. The alternating minimization algorithm and the  summable convergence are presented in Section <ref>. Simulation results on the IEEE 123-bus system and a real utility feeder are shown in Section <ref>. Section <ref> concludes the paper.





§ PROBLEM FORMULATION


This section first reviews the standard matrix completion formulation, and then  introduces how to build the data matrix using the time-series data in multiphase distribution systems. Next, the linearized power flow constraints are introduced. Finally,  the full matrix completion model for state estimation in multiphase distribution systems is derived, which combines the standard matrix completion model with the linear power flow constraints.


 §.§ Standard Matrix Completion





Let M∈ℝ^m× n  be the matrix we wish to reconstruct, and let Ω denote the set of indices for known elements of M. Define the observation operator P_Ω:ℝ^m× n→ℝ^m× n as:
   P_Ω(M)_ij=M_ij for (i,j) ∈Ω; 0 otherwise.

 The goal is to recover M from P_Ω(M),  under the assumption that M is low rank. Because  the rank function is nonconvex and the tightest convex  relaxation of  the  rank  function is the nuclear norm function  <cit.>, we consider the following regularized matrix completion model:
 If the low-rank structure of the matrix is true, we want to find the matrix with the minimum rank that is consistent with the observations, leading to:

    X∈ℝ^m× nminimize      rank(X) 
       subject to      P_Ω(X)=P_Ω(M),

where X is the decision variable and rank(X) is the rank of  X.

Unfortunately, this formulation is of no practical use because the rank minimization is NP-hard. One approach, proposed by <cit.>, is to  minimize the nuclear norm instead, which is the tightest  convex relaxation of the rank function. In addition, the measurements are susceptible to noise, and <cit.> proposed a noise-resilient constraint ||𝒫_Ω(X)-𝒫_Ω(M)||_F^2≤δ (δ>0).
Denote the nuclear norm of X by 

||X||_*=∑_i=1^min{m,n}σ_i(X),

    ||X||_*=∑_i=1^min{m,n}σ_i(X),

where σ_i(X) is the ith largest singular value of X.  
Then the new model with convex relaxation and noise-resilience is  formulated as:

    X∈ℝ^m× nmin       ||X||_*+μ/2||P_Ω(X)-P_Ω(M)||_F^2,

where X is the decision variable, ||X||_*=∑_i=1^min{m,n}σ_i(X) is the nuclear norm of X with σ_i(X) denoting the ith largest singular value of X, and μ>0 is a weight parameter.



 §.§ Matrix Formulation  in Multiphase Distribution Systems


For the scope of this paper, 

we consider a general three-phase distribution network consisting of one slack bus and a given number of multiphase PQ buses. 


Matrix completion was shown  effective for state estimation in general three-phase distribution systems by <cit.>, which formulated the data matrix M using  single-snapshot information. Considering  the temporal correlation of the data, we set up the data matrix M using time-series data. 



 Assume that the voltage phasor and other measurements at the slack bus are known. So we use only  measurements at nonslack buses to form the data matrix. Let 𝒫 denote the set of phases at all nonslack buses and |𝒫| be the number of all phases at  nonslack buses. 
 The measurements we  use in this matrix are the real and imaginary parts of voltage phasor, voltage magnitude, and the active and reactive power injection at each phase of nonslack buses.  Consider a time series t=1,…, T. Let M^t denote the measurement matrix at time t such that each column represents a phase and each row represents a quantity relevant to the phase. To be specific, for each phase i∈𝒫,
the corresponding column of M^t is of the form:

    [(v_i), (v_i), |v_i|, (s_i),  (s_i)]^⊺,

where (·) and (·) are the real part and imaginary part of a complex variable, respectively, ⊺ is the transpose notation, v=[v_1,…, v_|𝒫|]^⊺∈ℂ^|𝒫| is the vector containing voltage phasors at each phase of nonslack buses, and s=[s_1,…, s_|𝒫|]^⊺∈ℂ^|𝒫| is the vector of power injections at each phase of nonslack buses. 
Then the matrix M is constructed by 

    M=[M^1; M^2; ⋯;  M^T] ∈ℝ^m × n,

where  T is the length of the time  series, and m = 5T, n=|𝒫|.
Note that the data matrix is not limited to the structure proposed above; it can accommodate any available measurements if these measurements are correlated such that M has the (approximate) low-rank property.


Because the rows and/or columns of the matrix M are related by temporal correlation, spatial correlation, or power flow equations, such a matrix has low rank. Fig. <ref> shows that the low-rank property holds for the IEEE 123-bus  system with one-time step, two-time step, and three-time step  data formulation.


  


 §.§ Linear Power Flow Constraints
 

As in <cit.>,
we use the linear approximation  of  power  flows  as  constraints  to  characterize the correlation among different measurement types, which provides  more  information  about  the  system,  and  thus  likely  requires fewer  measurements  to  recover  the  matrix  than the standard  matrix  completion  formulation.  





Recall that we use  𝒫, v∈ℂ^|𝒫|, and s∈ℂ^|𝒫|  to denote the set of phases, the vector of voltage phasors, and the vector of power injections at  all nonslack buses, 
respectively. 
 We employ approximations of the form:
 

    v      ≈ B [ (s); (s) ] +   w,  
    
    |v|    ≈ C [ (s); (s) ]  +|w|,





where the coefficients B,C∈ℂ^|𝒫|× 2|𝒫|, w∈ℂ^|𝒫| are derived by the power flow linearization methods developed in, e.g., <cit.>.

Write (<ref>) in the following equivalent form with A_1,A_2,A_3,A_4,C_1,C_2∈ℝ^|𝒫|× |𝒫|:
 

    (v)      ≈ [A_1   A_2] [ (s); (s) ] +  (w),  
    (v)      ≈ [A_3   A_4] [ (s); (s) ] +  (w),  
    
    |v|    ≈ [C_1   C_2] [ (s); (s) ]  +|w|,


where [0.9][A_1    A_2] = (B),[A_3    A_4]=(B), and [0.9][C_1   C_2]=C.

 We use v^t, s^t to denote the corresponding v,  s in (<ref>) at time slot t. Assume that v_0 (the slack bus voltage phasor) and the system topology are the same at different time steps, and then by using the linearization method in  <cit.>, A_1, A_2, A_3, A_4, C_1, C_2 are the same at different time steps.  Hence, we have the linear approximation of voltage phasor and magnitude at time t in (<ref>) with v,s being replaced by v^t and s^t.
 as: don't need to repeat this equation, just refer to (4)
  

    (v^t)      ≈ [A_1   A_2] [ (s^t); (s^t) ] +  (w),  
    (v^t)      ≈ [A_3   A_4] [ (s^t); (s^t) ] +  (w),  
    
    |v^t|    ≈ [C_1   C_2] [ (s^t); (s^t) ]  +|w|.




 For simplicity of expression in the sequel, we use the following model to express this linear model at time t=1,…, T: 
 
    y≈ Ax+b,

where: 
    y = [ (v^1); (v^1); |v^1|;     ⋮; (v^T); (v^T); |v^T| ],  A = [ A_1 A_2   0   ⋯   0   0; A_3 A_4   0   ⋯   0   0; C_1 C_2   0   ⋯   0   0;   ⋮   ⋮   ⋮   ⋱   0   0;   0   0   0   ⋯ A_1 A_2;   0   0   0   ⋯ A_3 A_4;   0   0   0   ⋯ C_1 C_2 ],


    x= [ (s^1)^⊺ (s^1)^⊺ … (s^T)^⊺ (s^T)^⊺ ]^⊺,
 and
     
    b =[ (w)^⊺ (w)^⊺ |w|^⊺ … (w)^⊺ (w)^⊺ |w|^⊺ ]^⊺


    x= [ (s^1); (s^1);     ⋮; (s^T); (s^T) ],
     and b =
       [ (w); (w); |w|;   ⋮; (w); (w); |w| ]

with A∈ℝ^3T|𝒫|× 2T|𝒫|,x∈ℝ^2T|𝒫|, and y, b∈ℝ^3T|𝒫|.


The linear model (<ref>) will be incorporated into the matrix completion formulation in the following subsection.



 §.§ The Regularized Matrix Completion for State Estimation



By incorporating the  linear power flow model (<ref>) into  problem (<ref>) as a regularization term, we obtain the following optimization problem:

    X∈ℝ^m× nminimize       ||X||_*
       subject to       ||𝒫_Ω(X)-𝒫_Ω(M)||_F^2≤δ,
              |||y-(Ax+b)||_2^2≤β,

where β≥ 0 is a parameter that can be chosen based on the extent of approximation.

 
 
 I suggest the following:
 
     
  * Stick with problem (9) everywhere, if this is the problem we are solving.
     
  * State this problem already in Section II.A, with general regularization term.
     
  * I reformulated a bit  problem (9) below – check if this makes sense. We don't need to define X again, it was already defined above. 
 
 
 
    min    ||X||_*+μ/2||P_Ω(X-M)||_F^2 +ν/2||y-(Ax+b)||_2^2 
       over    X ∈ℝ^m× n,y ∈ℝ^3/5mn,x ∈ℝ^2/5mn,
       s.t.   [0.95]y = [
        a_1^⊺X    
        a_2^⊺X   
        a_3^⊺X … 
        a__3T-2^⊺X  
        a__3T-1^⊺X  
        a__3T^⊺X 
    ]^⊺, 
          [0.95]x = [
       c_1^⊺X  
         c_2^⊺X  … 
        c__2T-1^⊺X  
        c__2T^⊺X
             ]^⊺,
 
where, ν>0 is a weight parameter,  m = 5T, n = |𝒫|,
a_3(t-1)+i=e_5(t-1)+i and c_2(t-1)+i=e_5(t-1)+3+i are standard basis vectors in ℝ^m, in which t denotes the time step number and i denotes the ith a or c at the t-th time step. For example, if a_1, t =1, i =1, then, a_1 = e_1∈ℝ^m; if c_1, t=1, i = 1, then c_1=e_4∈ℝ^m. 

 
 
 We consider an approximate way to solve problem (<ref>)
by solving the following  unconstrained problem:

    X ,y,xminimize     ||X||_*+μ/2||𝒫_Ω(X)- 𝒫_Ω(M)||_F^2
        +ν/2||y-(Ax+b)||_2^2.

where    μ,ν>0 are the penalty parameters, y and x are variables mapped from X, and the specific mappings   are:

    X = [ (v^1)^⊺; (v^1)^⊺; |v^1|^⊺; (s^1)^⊺; (s^1)^⊺;       ⋮; (v^T)^⊺; (v^T)^⊺; |v^T|^⊺; (s^T)^⊺; (s^T)^⊺ ], y = [    X^⊺a_1;    X^⊺a_2;    X^⊺a_3;         ⋮; X^⊺a_3T-2; X^⊺a_3T-1;   X^⊺a_3T ],x = [    X^⊺c_1;    X^⊺c_2;         ⋮; X^⊺c_2T-1;   X^⊺c_2T ]

with X∈ℝ^m× n, y∈ℝ^3/5mn, and x∈ℝ^2/5mn.


    a_1    = [ 1 0 0 0 0 0 ⋯ 0 0 0 ]^T∈ℝ^5T,
    
        a_2    = [ 0 1 0 0 0 0 ⋯ 0 0 0 ]^T∈ℝ^5T,
    
        a_3    = [ 0 0 1 0 0 0 ⋯ 0 0 0 ]^T∈ℝ^5T,
        ⋮
    
        a_3T-2   = [ 0 0 0 0 0 0 ⋯ 1 0 0 ]^T∈ℝ^5T,
    
        a_3T-1   = [ 0 0 0 0 0 0 ⋯ 0 1 0 ]^T∈ℝ^5T,
    
        a_3T   = [ 0 0 0 0 0 0 ⋯ 0 0 1 ]^T∈ℝ^5T,
    
        c_1     =  [ 0 0 0 1 0 0 ⋯ 0 0 0 ]^T∈ℝ^5T,
    
         c_2     =  [ 0 0 0 0 1 0 ⋯ 0 0 0 ]^T∈ℝ^5T,
        ⋮
    
          c_2T-1    =  [ 0 0 0 0 0 0 ⋯ 0 1 0 ]^T∈ℝ^5T,
    
          c_2T    =  [ 0 0 0 0 0 0 ⋯ 0 0 1 ]^T∈ℝ^5T.











§ ALTERNATING MINIMIZATION FOR SOLVING THE FULL MATRIX COMPLETION MODEL


This section first  proposes the alternating minimization algorithm for the full matrix completion model in multiphase distribution system state estimation, and then proves the  summable convergence of the proposed algorithm.


 §.§ Alternating Minimization Algorithm

Any matrix X∈ℝ^m× n of a rank up to r has a matrix product X= UV form with U∈ℝ^m× r and V∈ℝ^r× n, and its nuclear norm can be expressed by the Frobenius norm of U and V as follows <cit.>:


    ||X||_*:=    U ∈ℝ^m × r,V ∈ℝ^r × nmin 1/2 (||U||_F^2+||V||_F^2) 
           s. t.        
     X=UV.


Note that for sufficiently small r (r≪min{m,n}), using this characterization of the nuclear norm  allows us to dramatically reduce the size of the problem.  We now use (<ref>) to reformulate (<ref>) as follows:

    [0.95]U ∈ℝ^m × r,
    V ∈ℝ^r × nmin   [0.95]1/2 (||U||_F^2+||V||_F^2)+μ/2 ||P_Ω(UV-M)||_F^2
                  [0.95]+ν/2||f_1(UV)-(Af_2(UV)+b)||_2^2,
       s.t.   [0.86]f_1(UV) = [
        a_1^⊺X    
        a_2^⊺X   
        a_3^⊺X  … 
        a__3T-2^⊺X  
        a__3T-1^⊺X  
        a__3T^⊺X 
    ]^⊺, 
          [0.86]f_2(UV) = [
       c_1^⊺X  
         c_2^⊺X  … 
        c__2T-1^⊺X  
        c__2T^⊺X
             ]^⊺,
          [0.86]X =UV.



The alternating minimization algorithm that updates the variables in an alternating fashion is proved to be one of the most accurate and effective methods for solving the matrix completion problem <cit.>. For our matrix completion formulation in state estimation, the alternating minimization algorithm solves for U and V alternatively while fixing the other factor.
The pseudo-code of the alternating minimization algorithm for the formulation (<ref>) is given in Algorithm <ref>.






 §.§ Convergence of the Alternating Minimization Algorithm




Because of the bilinear form X=UV, problem (<ref>) is not jointly convex with respect to any variable.
In general, there is no guarantee that a stationary point of a nonconvex problem will coincide with the global optimum. However, <cit.> proved the  summable convergence of the alternating minimization algorithm if the objective function satisfies the following four conditions: 1) it  is  continuous  in  its  domain  and  always bigger  than  negative  infinity;  2) it  has  a Nash equilibrium; 3) it is strongly convex with respect  to  each  variable; and 4) it satisfies the Kurdyka-Łojasiewicz property (please see Definition 2.5 in  <cit.> for the definition). We show below that our objective function satisfies these four conditions;  thus, we have the following proposition. 



The alternating minimization algorithm (Algorithm 1) satisfies the summable convergence, i.e.: 
    ∑_k=1^+∞||U^(k+1)-U^(k)||_F+||V^(k+1)-V^(k)||_F<+∞.
 



    ||U^(k+1)-U^(k)||_F^2+||V^(k+1)-V^(k)||_F^2≤
    
        ||f(U^(k+1),V^(k+1))-f(U^(k),V^(k))||^2


The objective function f can be written as follows:

    f(U,V) =    1/2 (||U||_F^2+||V||_F^2)+μ/2 ||P_Ω(UV-M)||_F^2
        +ν/2||f_1(UV)-(Af_2(UV)+b)||_2^2,

where:

    ||f_1(UV)-(Af_2(UV)+b)||_2^2
    
        =   [0.9]∑_i=1^T||X^⊺a__3i-2-(A_1X^⊺c__2i-1+A_2X^⊺c__2i)-(w)||_2^2
         [0.9]+∑_i=1^T||X^⊺a__3i-1-(A_3X^⊺c__2i-1+A_4X^⊺c__2i)-(w)||_2^2
         [0.9]+∑_i=1^T||X^⊺a__3i-(C_1X^⊺c__2i-1+C_2X^⊺c__2i)-|w|||_2^2

with X=UV.

1) The objective function f(U,V) is continuous in its domain, and inf_U,V∈domainf(U,V)>-∞.

2) Problem (<ref>) has a Nash point.

3) The objective function f(U,V) is strongly convex with respect to U and V, respectively.



First, the domain of the objective function is (𝒰,𝒱)={(U,V):U∈ℝ^m× r, V∈ℝ^n× r}.
By the form of f(U,V) in (<ref>), we have that f(U,V) is continuous in its domain (𝒰,𝒱), and it is non-negative for any (U,V)∈ (𝒰,𝒱) .



Second,  to prove that (<ref>) has a Nash point, we prove that (<ref>) has a global minimizer, which by definition is a Nash point.  By  the definition of f(U,V) in (<ref>), the following expression holds:

    lim_||U||→ +∞,
    
    ||V||→ +∞ f(U,V)=+∞,

which means that f(U,V) is coercive <cit.>.
 Then, by Theorem 2 in <cit.>, problem (<ref>) has a global minimizer.

Third, we prove that f(U,V) is strongly convex with respect to U and V, respectively. 
To prove that f(U,V) is strongly convex with respect to U, we prove that ∂^2f(U,V)/∂ U^2≻ 0 as follows.

Let g(U):=||f_1(UV)-(Af_2(UV)+b)||_2^2 for a fixed V. Then we have 


    ∂^2f(U,V)/∂ U^2= I_r× r⊗ I_m× m+μ V^TV⊗ I_m× m+ν/2∂^2g/∂ U^2.




    ∂^2f(U,V)/∂ U^2=    I_mr× mr+μ (V⊗ I_m× m)^T R_Ω (V⊗ I_m× m) 
           +ν/2∂^2g/∂ U^2,

where R_Ω∈ℝ^mr× mr is the diagonal matrix such that 
(R_Ω)_ii=1 if  ( ( i  n) +1, ⌈i/m⌉) ∈Ω and 0 otherwise.

By the form of g in (<ref>), we  focus only at the second-order partial derivative of the first term in the first summation term in (<ref>); treatment of the other terms can be derived similarly. Let g_1  denote the first term in the first summation term in (<ref>).
Then

    ∂^2g_1/∂ U^2   =
        ∂^2(c_1^TUV^TA_1^TA_1VU^Tc_1)/∂ U^2+∂^2(c_2^TUV^TA_2^TA_2VU^Tc_2)/∂ U^2
        +∂^2(a_1^TUV^TVU^Ta_1)/∂ U^2+∂^2(2c_1^TUV^TA_1^TA_2VU^Tc_2)/∂ U^2
          -∂^2(2c_1^TUV^TA_1^TVU^Ta_1)/∂ U^2-∂^2(2c_2^TUV^TA_2^TVU^Ta_1)/∂ U^2
       = 2V^TA_1^TA_1V⊗ c_1c_1^T+2V^TA_2^TA_2V⊗ c_2c_2^T
        +2V^TV⊗ c_1c_1^T+2V^TA_1^TA_2V⊗ c_1c_2^T
        +2V^TA_2^TA_1V⊗ c_2c_1^T-2V^TA_1^TV⊗ c_1a_1^T
        -2V^TA_1V⊗ a_1c_1^T-2V^TA_2^TV⊗ c_2a_1^T
        -2V^TA_2V⊗ a_1c_2^T.

It is hard to prove that ∂^2g_1/∂ U^2≽ 0_mr× mr by the form of ∂^2g_1/∂ U^2 in the above equation. 
We will prove that 
g_1 is a convex function by definition, then using that g_1:ℝ^m× r→ℝ is convex if and only if ∂^2g_1/∂ U^2≽ 0 to derive that  ∂^2g_1/∂ U^2≽ 0.

By (<ref>), we have that:

    g_1(U) = ||A_1V^⊺U^⊺c_1+A_2V^⊺U^⊺c_2-V^⊺U^⊺a_1+(w)||_2^2.

Because 𝒰 denotes the set of all matrices of size m× r, we have that 𝒰 is a convex set.
Then by the form of g_1(U), it is easy to derive that  ∀ U_1, U_2∈𝒰, ∀ t∈ [0,1]: 
    g_1(tU_1+(1-t)U_2)≤ tg_1(U_1)+(1-t)g_1(U_2).


    g_1(tU_1+(1-t)U_2)
       =||tA_1V^⊺U_1^⊺c_1+(1-t)A_1V^⊺U_2^⊺c_1+tA_2V^⊺U_1^⊺c_2
          +(1-t)A_2V^⊺U_2^⊺c_2-tV^⊺U_1^⊺a_1-(1-t)V^⊺U_2^⊺a_1
          +t(w)+ (1-t)(w)||_2^2
       ≤ t||A_1V^⊺U_1^⊺c_1+A_2V^⊺U_1^⊺c_2-V^⊺U_1^⊺a_1+(w)||_2^2
        [0.95]+ (1-t)||A_1V^⊺U_2^⊺c_1+A_2V^⊺U_2^⊺c_2-V^⊺U_2^⊺a_1+(w)||_2^2
       =tg_1(U_1)+(1-t)g_1(U_2).






Hence, g_1 is convex, implying that ∂^2g_1/∂ U^2≽ 0.

Because all the other terms in (<ref>) share the similar form with g_1, we can similarly prove that all the other terms in (<ref>) are convex functions with respect to U; thus, we have that ∂^2g/∂ U^2≽ 0.
By (<ref>), we have that 
I_mr× mr≻ 0, and μ (V⊗ I_m× m)^T R_Ω (V⊗ I_m× m) ≽ 0; therefore, for μ>0, ν>0, ∂^2g/∂ U^2≽ 0, we have
∂^2f(U,V)/∂ U^2≻ 0.

Hence, we conclude that f is strongly convex with respect to U. The proof that f is strongly convex with respect to V is similar and is omitted for brevity.


By <cit.>, any real analytic function satisfies the Kurdyka-Łojasiewicz  property. And by definition of the real analytic function, the objective function (<ref>) is real analytic; thus, (<ref>) satisfies the Kurdyka-Łojasiewicz property.

This concludes the proof.


We can actually do better than square summable.  this section proves that our objective meets the assumptions used in Yangyang Xu's paper, which utilizes the KL property to further show that the norm difference of the iterates is summable 
April, that will be great if we can have stronger convergence results. But what I am concerned is whether you can write down the proof this week. Because I am working on the simulation results and we'd better finish the paper writing by the end of this week (the submission deadline is 9/23).
You've already proved it.  The paper you cite proves iterates are summable in addition to being square summable, so you can just change the wording in the last sentence


§ SIMULATION RESULTS


This section demonstrates Algorithm <ref> on the IEEE 123-bus test system <cit.> and  a real utility feeder system. The two systems are both three-phase unbalanced radial distribution systems, in which buses are one-, two-, or three-phase. The slack bus is three-phase, and the total number of phases at all buses for the IEEE 123-bus system is 263 and that for the real utility feeder system is 1234. The data for both systems were simulated at 1-minute resolution using power flow analysis with diversified load and solar profiles that were created for each bus using realistic solar irradiance and load consumption data.
The voltage magnitudes range from 0.95 to 1 p.u., and the voltage angles are around 0 or ± 120 degrees.



As in Section <ref>, the data matrix is formulated using  the  real and imaginary parts of voltage phasor,  voltage  magnitude,  and active  and reactive power injection at each phase of nonslack buses. We consider  the real and imaginary parts of voltage phasor as variables, and the voltage magnitude,  active power, and reactive power as potentially known measurements. In our experiments, we define the percentage of known measurements as the ratio of the number of known measurements and the total number of voltage magnitude, active power, and reactive power in the data matrix. If less than two thirds of the measurements are known, there are fewer measurements than variables, under which case, the system is underdetermined and the WLS cannot generate a unique solution. For all of our experiments,  the known measurements are added with noise with zero mean and 1% of the true value as the  standard deviation. And the simulation results are based on the average of 50 runs for each scenario.


 §.§ Performance on IEEE 123-Bus System

First, we implement our algorithm with one-time step data formulation at one time slot when solar generations are nonzero with 10% to 70% measurements available (where WLS is under-determined). We compare
the alternating minimization algorithm with the SDP solver used by <cit.>. Both algorithms are run with MATLAB on a laptop with 1.9-GHz CPU and 32 GB RAM. 

 Fig. <ref> shows one representative convergence result (with 50% available measurements)  of the alternating minimization algorithm. 
Fig. <ref> shows that the mean absolute percentage errors (MAPEs) of the voltage magnitude  and the mean absolute errors (MAEs) of the voltage angle are comparable for the  alternating minimization algorithm and the SDP solver. The running time for the  alternating minimization algorithm is about 15 seconds  for different percentages of available measurements, whereas that for the SDP solver is more than 70 seconds. These results show that the  alternating minimization algorithm estimates voltage phasors more  computationally efficient than the SDP solver under similar performance. The running time for the SDP solver increases  as more measurements are available. This is  because the SDP solver solves a constrained SDP problem, where more measurements being available means more constraints to be satisfied, and thus longer time is required to achieve the convergence.





To test the computation  time  and  estimation  accuracy as more data are used, we  implement the  alternating minimization algorithm on the two-time step and three-time step data matrix formulation with 10%, 30%, 50%, and  70% measurements available. Fig. <ref> shows that the MAPEs of the voltage magnitudes and the MAEs of the voltage angles for the two-time step and three-time step cases are much better than the one-time step case. This is because the matrix completion problem minimizes the rank, and the larger the time steps are, the smaller the rank of the matrix is compared to its size (which is shown in Fig. <ref>) and the more accurate the estimates are. The  running time also increases with the length of time steps, however, because of the increased size of the matrices.
The MAPEs of the voltage magnitudes and MAEs of the voltage angles for 
the three-time step case is slightly better than the two-time step case overall, but the running time almost doubles. However, the running time for the three scenarios shows that the alternating minimization algorithm for state estimation can be implemented in real time for the IEEE 123-bus system.





 §.§ Performance on A Real Feeder System

In order to demonstrate the scalability of our algorithm, we implement it on a real utility feeder system with 1234 phases under the assumption that 50%  measurements are available  for one-time, two-time, and three-time step matrix formulation, respectively. The performance in terms of MAPE, MAE, and running time are shown in Table <ref>. The results show that the performance improves with more data being used for building the matrix and the computation time also increases. However, the running time shows that the algorithm can be applicable in real-time estimation.



§ CONCLUSIONS


This paper applied the constrained matrix completion formulation for state estimation in multiphase distribution systems. Different from the constrained matrix completion, the data matrix was formulated using time-series data and  the alternating minimization algorithm was applied to solve the  constrained matrix completion. The  summable convergence of the alternating minimization algorithm was proved. In addition, the efficacy and scalability of the algorithm were demonstrated via the IEEE 123-bus system and a real utility feeder system.  
 
 
 

IEEEtran


