Unified value-based feedback, optimization and risk management in complex electric energy systems

Abstract

The ideas in this paper are motivated by an increased need for systematic data-enabled resource management of large-scale electric energy systems. The basic control objective is to manage uncertain disturbances, power imbalances in particular, by optimizing available power resources. To that end, we start with a centralized optimal control problem formulation of system-level performance objective subject to complex interconnection constraints and constraints representing highly heterogeneous internal dynamics of system components. To manage spatial complexity, an inherent multi-layered structure is utilized by modeling interconnection constraints in terms of unified power variables and their dynamics. Similarly, the internal dynamics of components and sub-systems (modules), including their primary automated feedback control, is modeled so that their input–output characterization is also expressed in terms of power variables. This representation is shown to be key to managing the multi-spatial complexity of the problem. In this unifying energy/power state space, the system constraints are all fundamentally convex, resulting in the convex dynamic optimization problem, for typically utilized quadratic cost functions. Based on this, an interactive multi-layered modeling and control method is introduced. While the approach is fundamentally based on the primal–dual decomposition of the centralized problem, this is formulated for the first time for the couple real-reactive power problem. It is also is proposed for the first time to utilize sensitivity functions of distributed agents for solving the primal distributed problem. Iterative communication complexity typically required for convergence of point-wise information exchange is replaced by the embedded distributed optimization by the modules when creating these functions. A theoretical proof of the convergence claim is given. Notably, the inherent multi-temporal complexity is managed by performing model predictive control (MPC)-based decision making when solving distributed primal problems. The formulation enables distributed decision-makers to value uncertainties and related risks according to their preferences. Ultimately, the distributed decision making results in creating a bid function to be used at the coordinating market-clearing level. The optimization approach in this paper provides a theoretical foundation for next-generation Supervisory Control and Data Acquisition (SCADA) in support of a Dynamic Monitoring and Decision Systems (DyMonDS) for a multi-layered interactive market implementation in which the grid users follow their sub-objectives and the higher layers coordinate interconnected sub-systems and the high-level system objectives. This forms a theoretically sound basis for designing IT-enabled protocols for secure operations, planning, and markets.

Appendices

Appendix

1.1 A Theoretical basis for the provability of DyMonDS-based bids

Let us first pose the problem for allocation of resources in general form. From the benchmark problem described in Eq. (19), one can see that the coupling constraints are all linear in interaction space. Let us use the notation $z_i[k]$ to represent interactions including $E_i, p_i, P_i, \dot{Q}_i$ for all time samples $kT_t \in [t, t+H]$. The internal variables of including $x_i[k], u_i[k], m_i[k]$ for the same time samples are all abstracted in the vector $x_i$ (with slight abuse of notation). Let x, z denote the vectors $x_i, z_i$ for all components stacked up.

The centralized problem can then be posed in a general form as shown below:

$$\begin{aligned} \begin{array}{l} {(C):}\,\,{\mathop {\min }\limits _{{z_i}}} \sum \limits _{i \in \mathcal {N}} {f_i(z_i)} \\ s.t.\\ \quad\quad g(z) \le 0 \quad \left( \lambda \right) \\ \qquad h_i(x_i,z_i) \le 0 \quad \left( \mu _i\right) \\ \end{array} \end{aligned}$$

(36)

The first constraint function g(z) is a coupling constraint while the second one $h_i(x_i,z_i)$ is specific to each agent i. Solving the problem in Eq. (36) turns out to be complicated since the number of components given by $|\mathcal {N}|$ is extremely large. While the coupling constraints given by g(z) are linear, the number of internal constraints modeled through $h_i$ for each agent i can be high, leading to a large-scale problem to be solved at once.

In this general formulation and the proofs to follow, a few assumptions are made:

Assumption 1

The cost function $f_i: \mathbb {R}^{n_{i}} \rightarrow \mathbb {R}$ is a smooth convex function which is continuously differentiable with Lipschitz continuous gradients $L(f_i) > 0$ given as

$$\begin{aligned} \left\| \nabla {f_i}({x_i}) - \nabla {f_i}({y_i})\right\| \le L(f_i)\left\| {x_i} - {y_i}\right\| \quad \forall x_i,y_i\in \mathbb {R}^{n_i} \end{aligned}$$

(37)

where $\Vert .\Vert$ denotes the standard Euclidean norm.

Assumption 2

The constraint set given in (36) is a closed convex set.

Exploiting the linearity of the coupling constraints, we have proposed a distributed approach towards solving the problem (C). As explained in Sect. 6, we propose to have the system operator solve the master problem comprising the coupling constraints alone which are of the same order as the number of nodes in the network. The increased cardinality of the problem created by local constraints and their preferences are rather abstracted through bid functions as computed by solving the agent-level problem in Eq. (38)

$$\begin{aligned} \begin{array}{l} {(Agent-i)}:{\mathop {\min }\limits _{{x_i},z_i}} {C_i}({z_i}) + {\lambda }_j^T g_j(z_i)\\ \qquad \qquad s.t.\\ \qquad \qquad h_i(x_i, z_i) \le 0 \\ \end{array} \end{aligned}$$

(38)

Here, $\lambda _j$ corresponds to the vector of all the elements of vector $\lambda$ for which there exists dependence of the constraint on the elements of $z_i$ i.e. $j \in \mathcal {J}$ to denote constraint indexes to which $z_i$ contributes.

By solving the problem $\left( Agent-i\right)$ for perturbations in each of the elements in $\lambda _j$, one can obtain the bid functions $B^j_i(z_i)$ for use by the operator to solve the system level problem. The willingness of the agent to participate depends on its internal constraints but also on the price incentive. The bid is such that it makes marginally the same amount of cost as it gets paid through the respective $\lambda _j$. As a result,

$$\begin{aligned} \frac{\partial B_i^j}{\partial z_i} = -\lambda _j = a^j_iz_i + b^j_i \end{aligned}$$

(39)

The coefficients $a^j_i$ and $b^j_i$ are computed empirically by obtaining the points $z_i$ by solving $\left( Agent-i\right)$ for slight perturbations in $\lambda$. This function is integrated with respect to $z_i$ to obtain

$$\begin{aligned} B_i(z_i) = \sum _{j \in \mathcal {J}} B_i^j(z_i) = \frac{1}{2}z_i^T a_i z_i + b_i z_i \end{aligned}$$

(40)

Here, $a_i = \sum _{j \in \mathcal {J}}a_i^j$ and $b_i = \sum _{j \in \mathcal {J}}b_i^j$. Along with the bid functions, the feasibility region of the agent i is projected onto the $z_i$ plane represented by the space $\varPi _{i}$. This projected space varies with time as operating conditions change. This plane is found by utilizing the constraints Eq. (19f). These bid functions along with the time varying limits of the feasible region of interaction variables are all collected by the system operator to solve the following problem:

$$\begin{aligned} \begin{array}{*{20}{l}} (S): {\mathop {\min }\limits _{z_i}}\sum \limits _{i \in \mathcal {N} } {{B_i}(z_i)}\\ {s.t.}\\ \qquad \qquad g(z) \le 0 \qquad \left( \lambda \right) \\ \qquad \qquad z_i \in \varPi _{g,i} \end{array} \end{aligned}$$

(41)

Assumption 3

The problem $\left( C\right)$ is solvable and there exists a unique minimizer $x^*$.

Notice that the projected space $\varPi _{i}$ is also convex since the projection map is a linear operator, the constraint set of problem (S) is convex. Since the bid functions are by construction convex, there would be a unique minimizer of $\left( S\right)$. We next prove that the optimal point obtained by the proposed distributed scheme is the same as the one obtained by the centralized approach.

1.2 A.1 Proof for approaching the centralized solution optimum operating point through proposed primal–dual decomposition strategy

Theorem 1

Under assumptions 1–2, the solution obtained by solving$\left( S\right)$results in the same solution as that of the$\left( C\right)$when the bids$B_i(z_i)$are created by solving problem in$\left( Agent-i\right)$in response to small perturbations around$\lambda = \lambda ^*$where$\lambda ^*$is the dual solution of the problem$\left( C\right)$.

Proof

Writing the KKT conditions for the problem $\left( C\right)$, we obtain:

Stationarity:

$$\begin{aligned}&\nabla _{z_i} C_i(z_i) + \lambda ^T \nabla _{z_i} g(z) + \mu _i^T \nabla _{z_i} h_{i} (z_i,x_i) = 0 \qquad \forall i \in \mathcal {N} \end{aligned}$$

(42a)

$$\begin{aligned}&\mu _i^T \nabla _{x_i} h_{i} (z_i,x_i) = 0 \end{aligned}$$

(42b)

Primal feasibility:

$$\begin{aligned} \qquad g(z)\le & {} 0 \end{aligned}$$

(42c)

$$\begin{aligned} h_{i}(x_i,z_i)\le & {} 0 \end{aligned}$$

(42d)

Dual feasibility:

$$\begin{aligned} \lambda\ge & {} 0 \qquad \,\,\, \end{aligned}$$

(42e)

$$\begin{aligned} \mu _{i}\ge & {} 0 \,\,\,\forall i\in \mathcal {N} \end{aligned}$$

(42f)

Complementarity slackness:

$$\begin{aligned} \lambda ^T g(z)= & {} 0 \end{aligned}$$

(42g)

$$\begin{aligned} \mu _i^T h_{i} (x_i,z_i)= & {} 0 \,\,\,\forall i \in \mathcal {N} \end{aligned}$$

(42h)

Let the solution of above set of equations be $x_i^*, z_i^*, \lambda ^*, \mu _i^*$.

Now, solving the problem posed in $(Agent-i)$ for a given value of $\lambda$, we obtain the solution of the agent-level problem. As explained previously, the bid function is constructed by integrating the empirically constructed sensitivity given by

$$\begin{aligned} \lambda (z_i) = -\nabla _{z_i} B_i(z_i) = a_i z_i + b_i \end{aligned}$$

(43)

This results in the bid function $B_i(z_i) = \frac{1}{2}a_iz_i^2 + b_iz_i$.

Analytically, one can write the KKT conditions of the problem $\left( Agent-i\right)$ to notice that the gradient of the bid function satisfies the following relation:

$$\begin{aligned} -\nabla _{z_i} B_i(z_i) \nabla _{z_i} g(z) + \nabla _{z_i} C_{i} (z_i) + \lambda ^T \nabla _{z_i}h(x_i,z_i) = 0 \end{aligned}$$

(44a)

while also satisfying the device-specific KKT conditions given by the set $S_i$ i.e.

$$\begin{aligned} S_i = \left\{ \left( x_i, z_i, \lambda _i, \mu _i\right) | \hbox {(42b)}, \hbox {(42d)}; \hbox {(42f)}; \hbox {(42h)} \right\} \end{aligned}$$

(44b)

The set above which when projected onto the space of the $z_i$ plane results in $z_i \in \varPi _{g,i}$, representing the space of projection of the half space represented by $g_i(x_i,z_i)$ onto the $z_i$ plane in order to commit operating conditions dependent limits of $z_i$ as a part of the bid function. Clearly, $z_i \in \varPi _{i} \implies \left( x_i, z_i, \lambda _i, \mu _i \right) \in S_i$ for some values of $x_i, \lambda _i, \mu _i$.

Now considering these bid functions and the projected constraint on $z_i$ in the problem $\left( S\right)$. The KKT conditions of $\left( S\right)$can be written as:

Stationarity:

$$\begin{aligned} \nabla B_i(z_i) + \lambda ^T \nabla _{z_i} g(z_i) \quad \forall i \end{aligned}$$

(45a)

Primal, dual feasibility and complementary slackness:

Equations (42c), (42e), (42g)

An added constraint of the bid function that $z_i$ belongs to the projections of the set $S_i$ ensures the other internal KKT constraints in Eq. (36) remain feasible. Upon substituting the relations in Eq. (44) into the stationarity condition in Eq. (45a) for each agent i and summing them up, one could obtain the KKT conditions in Eq. (42). Thus proved. $\square$

Remark 14

When the internal constraints $g_i, h_i$ are non-convex, multiple agent-level solutions may exist and thus centralized and distributed solutions may not be the same.

B Proof for the rate of convergence of the DyMonDS-based sensitivity bids

1.1 B.1 Agent-level bid function

Following the definitions and proof sketch in Beck and Teboulle (2009), we now analyse the convergence rates for the proposed approach. First, we define the indicator function for the constraint set of the problem $\left( Agent-i\right)$ as $H_i(x_i,z_i)$. With this, we obtain a function to be minimized at the agent-level as

$$\begin{aligned} F_i(x_i, z_i, \lambda _j) = f_i(x_i, z_i) + H_i(x_i,z_i) + \lambda _j^T g_j(z_i) \end{aligned}$$

(46)

Here, it is assumed that each agent’s cost function is strongly convex, smooth and differentiable with a Lipschitz constant $L_i(f_i)$.

Definition 1

The quadratic approximation of $F_i(x_i,z_i,\lambda )$ around the point $y_i$ for $L_i \ge L_i(f_i)$ is given as

$$\begin{aligned} Q_{L,i}(x_i,z_i, \lambda _j)= & {} f_i(x_i,y_i,\lambda _j) + \left\langle z_i-y_i, \nabla _{z_i} f_i(x_i,y_i,\lambda )\right\rangle \nonumber \\&+ \frac{L_i}{2} \Vert z_i -y_i\Vert _2^2 + H_i(x_i,z_i) + \lambda _j^Tg_j(z_i) \end{aligned}$$

(47)

It admits a unique minimizer

$$\begin{aligned} {p_{L,i}}(x_i,{y_i},\lambda _j ) = \mathop {\arg \min }\limits _{z_i} {Q_{L,i}}({x_i},{z_i},\lambda _j) \end{aligned}$$

(48)

By denoting $z_i$ as the minimizer and by using the notation $\gamma _i\left( x_i,y_i,\right) \in \frac{\partial H_i}{\partial z_i}(x_i,z_i)$, the optimality condition can be obtained from Eq. (48) by taking the partial derivative w.r.t $z_i$ and equating to zero.

$$\begin{aligned} \nabla _{z_i} f_i(x_i,y_i) + L_i\left( z_i-y_i\right) + \gamma _i\left( x_i,y_i\right) + \lambda _j^T\frac{\partial g_j}{\partial z_i} = 0 \end{aligned}$$

(49)

Since the coupling constraints are linear, the partial derivative $\frac{\partial g_j}{\partial z_i}$ is a constant matrix. As explained previously, the marginal bid function at the agent level is equal to the expression for $\lambda _j$, for which the analytical expression can be decomposed into a slope $a_i$ and intercept value $b_i$

$$\begin{aligned} \frac{\partial B^j_i}{\partial z_i} = -\lambda _j = \underbrace{\left( \frac{\partial g_j}{\partial z_i}\right) ^{\dagger }L_i}_{a^j_{i}} z_i + \underbrace{\left( \frac{\partial g_j}{\partial z_i}\right) ^{\dagger }\left( \nabla _z f_i(x_i,y_i) + \gamma _i\left( x_i,y_i\right) - Ly_i\right) }_{b^j_{i}} \end{aligned}$$

(50)

Here, $(.)^{\dagger }$ represents the pseudo inverse operation. The bid function then can be re-expressed as

$$\begin{aligned} B_i(z_i) = \sum _{j \in \mathcal {J}} B_i^j(z_i) = \frac{1}{2}z_i^T a_i(y_i) z_i + b_i(y_i) z_i \end{aligned}$$

(51)

By collecting these bid functions from all the agents, the system cost function is constructed as

$$\begin{aligned} f^s(z) = \sum _{i \in \mathcal {N}} = \frac{1}{2} z_i^Ta_i(y_i)z_i + b_i(y_i)^Tz_i \end{aligned}$$

(52)

Remark 15

Notice that the system-lvel cost function is convex by construction. At each of the agents, the respective quadratic functions constructed, approximate the variation of their cost and the internal constraints around y.

1.2 B.2 System-level optimality conditions

Similarly, defining the function $F^s\left( z\right)$ for the problem to be solved by the system operator as

$$\begin{aligned} F^s(z^k) = f^s(z^k,z^{k-1}) + G(z^k) \end{aligned}$$

(53)

where $G(z^k)$ represents the indicator function of the coupling constraints $g(z^k) \le 0$. We have $Q^s\left( z^k,y\right)$ to denote the quadratic approximation $F^s(z^k)$ around a point y for some $L^s>0$ Around an arbitrary point y, the cost function coefficients are computed using the initial guess $z^0$.

$$\begin{aligned} Q^s(z^k,y^k) = f^s(z^k,z^{k-1}) + \left\langle z^k-y, \nabla _z f^s\left( y,z^0\right) \right\rangle + \frac{L^s}{2} \Vert z^k-y^k\Vert ^2 + G(z^k) \end{aligned}$$

(54)

The unique minimizer of the quadratic approximation of $F^s(z^k)$ around a point y is denoted as $p_L^s\left( y\right)$. The optimality conditions can explicitly be written as

$$\begin{aligned} \nabla _z f^s(y,z^{0})+ L^s\left( z^k-y\right) + \gamma ^s\left( y\right) = 0 \end{aligned}$$

(55)

Let us now define the sequence $\left\{ z^k\right\}$ given by $z^k = p^s_L(z^{k-1})$

Remark 16

The series $F^s\left( z^k\right)$ is contracting.

Notice that $F^s(z^k) \le Q^s\left( z^k,z^{k-1}\right) \le Q^s\left( z^{k-1},z^{k-1}\right) = F^s(z^{k-1})$

We will now derive its convergence rate.

1.3 B.3 Convergence proof

We begin by revising one of the lemmas in Beck and Teboulle (2009) for the distributed setting here, which will be utilized in the rest of our convergence proof:

Lemma 1

For some$L^s>0$and$y\in \mathbb {R}^{n_z}$if$F^s\left( p_L^s\left( y\right) \right) \le Q_L^s\left( p_L^s\left( y\right) ,y\right)$. Then, $\forall z \in \mathbb {R}^{n_z}$,

$$\begin{aligned} F^s(z^k) - F^s\left( p_L^s\left( y\right) \right) \ge \frac{L^s}{2}\Vert p^s_L(y) - y\Vert ^2 + L^s\left\langle y-z^k, p^s_{L}(y)- y\right\rangle \end{aligned}$$

(56)

Proof

Notice that the condition in Lemma holds true if $L^s = L^s_y = diag(a_i(y_i))$ i.e. the matrix created by the quadratic coefficient of the bid function. Because of the convexity of $f^s$ with respect to $z^k$ and since G is an indicator function of the convex set of coupling constraints, we can establish the following relations:

$$\begin{aligned} f^s(z^k,z^{k-1})\ge & {} f^s(y, z^0) + \left\langle z^k - y, \nabla f^s(y, z^{0})\right\rangle \end{aligned}$$

(57)

$$\begin{aligned} g(z^k)\ge & {} g(p^s_L(y)) + \left\langle z^k - p^s_L(y), \gamma (y) \right\rangle \end{aligned}$$

(58)

We can now write the expression for the objective function as

$$\begin{aligned} F^s(z^k) \ge f^s(y, z^{0}) + \left\langle z^k - y, \nabla f^s(y, z^{0})\right\rangle + g(p_L(y)) + \left\langle z^k - p_L(y), \gamma (y) \right\rangle \end{aligned}$$

(59)

Furthermore, from the definition of Q, we have

$$\begin{aligned} Q(p_L(y),y) = f^s(y,y^{0}) + \left\langle p_L(y) - y, \nabla f^s(y, z^{0}) \right\rangle + \frac{L^s}{2} \Vert p_L(y) - y\Vert ^2+ g(p_L(y)) \end{aligned}$$

(60)

Subtracting the two equations, we obtain

$$\begin{aligned} F^s(z^k) - Q(p_L(y),y) \ge \left\langle z^k - p_L(y^k), \nabla f^s(y, z^{0}) + \gamma (y) \right\rangle - \frac{L^s}{2}\Vert p_L(y) - y\Vert ^2 \end{aligned}$$

(61)

Utilizing the assumption stated in the lemma and substituting the optimality condition in Eq. (55) in the right hand side, we obtain

$$\begin{aligned} F^s(z^k) - F^s(p_L(y))\ge & {} L \left\langle z^k - p_L(y^k), z^k - y \right\rangle - \frac{L^s}{2}\Vert p_L(y) - y\Vert ^2\nonumber \\= & {} \frac{L^s}{2}\Vert p_L(y) - y\Vert ^2 + L^s \left\langle y - z^k, p_L(y)-y \right\rangle \end{aligned}$$

(62)

$\square$

Theorem 2

Let$\left\{ z^k\right\}$denote the series produced by the update rule$z^k = p_L^s\left( z^{k-1}\right)$. Then, the system-level objective function converges at the rate of$\mathcal {O}\left( 1/k\right)$as follows:

$$\begin{aligned} F^s\left( z^k\right) -F^s\left( z^*\right) \le \frac{{\Vert z^0-z^*\Vert }_{L0}}{2k} \end{aligned}$$

(63)

where$z^*$is the optimal point and$L_0$is the diagonal matrix consisting of the quadratic coefficients of the bid functions of each of the agents.

Proof

Invoking the Lemma 1 with $z^k = z^*, y = z^n, L^s = L_n = diag(a_i(z^n_i))$, we have

$$\begin{aligned} {F^s}({z^*}) - {F^s}\left( {{z^{n + 1}}} \right)\ge & {} \frac{{{L_n}}}{2}\left( \Vert {{z^{n + 1}} - {z^n}\Vert ^2 + 2{L_n}\left\langle {{z^n} - {z^*},{z^{n + 1}} - {z^n}} \right\rangle } \right) \nonumber \\= & {} \frac{L_n}{2} \left( {\Vert {z^*} - {z^{n + 1}}\Vert }^2- {\Vert {z^*} - {z^{n}\Vert }}^2\right) \nonumber \\= & {} \frac{1}{2} \left( {\Vert {z^*} - {z^{n + 1}}\Vert ^2_{L_n}}- {\Vert {z^*} - {z^{n}\Vert ^2_{L_n}}}\right) \end{aligned}$$

(64)

Next, invoking the Lemma 1 with $z^k = z^n, y = z^n, L^s = L_n = diag(a_i(z^n_i))$, we have

$$\begin{aligned} {F^s}({z^n}) - {F^s}\left( {{z^{n + 1}}} \right) \ge \frac{1}{2} {\Vert {z^{n + 1}} - {z^n}\Vert ^2}_{L_n} \end{aligned}$$

(65)

Summing the inequality in Eq. (64) over $n = 0,1,\ldots k-1$, we have

$$\begin{aligned} k{F^s}({z^*}) - \sum \limits _{n = 0}^{k - 1} {{F^s}\left( {{z^{n + 1}}} \right) } \ge \Vert z^* - z^k\Vert ^2_{L_{k-1}} - \Vert z^* - z^0\Vert ^2_{L_{0}} \end{aligned}$$

(66)

Next, multiplying the inequality in Eq. (65) by n and summing the result over $n = 0,1,, \ldots k-1$, we have

$$\begin{aligned} -k{F^s}({z^k}) + \sum \limits _{n = 0}^{k - 1} {{F^s}\left( {{z^{n + 1}}} \right) } \ge \sum \limits _{n = 0}^{k - 1} {n \Vert z^{n+1} - z^0\Vert _{L_n}^2} \end{aligned}$$

(67)

Now, by adding equations (66) and (67), we obtain

$$\begin{aligned} F^s\left( z^k\right) -F^s\left( z^*\right) \le \frac{{\Vert z^*-z^0\Vert }_{L0}}{2k} \end{aligned}$$

(68)

$\square$

Here, the cleared values at each iteration are communicated to the agents which are utilized to compute the new bid functions. As a result, at $k^{th}$ iteration, the bids are computed using the solution obtain by system coordinator at $(k-1)^{th}$ iteration. The system level objective function convergence stated above indicates that these bid functions converge at a rate of $\mathcal {O}(1/k)$.

Note also that error at each iteration is because of the quadratic approximation of the bid function by the agents. The error in the This result combined with the existence of unique solution as obtained by the centralized and distributed solution strategies as shown in Theorem 1 ensures that the DyMonDS approach leads to the centralized optimal solution.

Remark 17

The computational complexity through this approach of minimal coordination is of the order $\mathcal {O}\left( n_z^3\right) + \sum _{i=1}^{|\mathcal {N}|} \mathcal {O}\left( n_{x_i}^3\right)$, where $n_z, n_{x,i}$ respectively represent the number of interaction variables, and number of state variables in each of the components of the network. Typically $n_{x_i}$ for each i is much smaller than $n_{z}$ of the entire system.

The centralized solution complexity would be of the order $\mathcal {O}\left( \left( n_z + \sum _{i=1}^{|\mathcal {N}|} n_{x_i}\right) ^3\right)$!

In order to compare with iterative approaches such as ADMM, simple update equations are typically utilized in addition to solving the local agent problems with a computational cost of $\sum _{i=1}^{|\mathcal {N}|}\mathcal {O}\left( n_{x_i}^3\right)$. However, the number of iterations needed is typically of the order $1e^6$, for obtaining an accuracy of $1e^{-3}$ in the solution vector.

1.4 B.4 Proof-of-concept numerical simulations

In order to show the effectiveness of the proposed DyMonDS-based interactive scheme to solving an optimization problem coupled across multiple agents, we consider a communication network shown in Fig. 8.

Fig. 8

Test network: Source destination pair is displayed with same color; $L_l: x_i, x_j \ldots$ denotes that the flow on link l is due to the flow from sources $i, j, \ldots$ (Beck et al. 2014)

It has two source-load pairs and the objective is to maximize the network utility. Each link has dedicated source power injections as shown in the Fig. 8 with an upper limit on the capacity of the flow through the wires. Let each of the sources have the utility functions $u_i(x_1) = C_ilog(x_1 + 0.1)$ where $C_1 = 10$ and $C_2 = 20$. Let the set of sources be denoted using the set $\mathcal {S}$ and let the set of sources utilizing the link l be denoted using $\mathcal {S}(l)$. Similarly, let the length of the path in use by source i be denoted as $\mathcal {L}(i)$. Letting the capacity of all the links be equal to 1, the network utility maximization (NUM) problem can be posed as

$$\begin{aligned} \begin{array}{l} {\mathop {\max }\limits _{x}} \sum \limits _{i \in \mathcal{S}} {{C_i}\log ({x_1} + 0.1)} \\ \sum \limits _{i \in \mathcal{S}(l)} {{x_i} \le 1\quad \forall l \in \mathcal{L}} \\ {x_i} \ge 0\quad \forall i \in \mathcal{S} \end{array} \end{aligned}$$

(69)

In order to apply DyMonDS-based method, this formulation is a degenerate case of the one posed in Eq. (36) where there are no internal constraints and all the variables appear the intersections. Thus, $x_i$ in this formulation is to be treated like $z_i$ in the formulation in Then the agent-specific formulation is to just optimize its own utility function given the lagrange multipliers corresponding to its coupling constraints. This can be posed as follows:

$$\begin{aligned} \begin{array}{l} {\mathop {\max }\limits _{x_i}} {{C_i}\log ({x_i} + 0.1)} - \sum _{j\in \mathcal {L}(i)} \lambda _j x_i \\ {x_i} \ge 0 \end{array} \end{aligned}$$

(70)

This problem is solved to obtain bids using the relations in Eq. (50) which reduces to the following for this degenerate case of having no internal constraints

$$\begin{aligned} \frac{\partial B_i}{\partial x_i} = \underbrace{\nabla ^2_xf_i(y)}_{a^{k}_i} x_i + \underbrace{\left( \nabla _x f_i(y) - \nabla ^2_x f_i(y) y\right) }_{b^k_i} \end{aligned}$$

(71)

Here $f_i(x_i) = C_i log(x_i + 0.1)$, y is the solution obtained by solving the agent level problem from $\lambda ^{k-1}$ in the previous iterate utilized to obtain the bid function for the next iteration. The bid constructed as $B^k_i(x_i) = \frac{1}{2}a_i^k x_i^2 + b_i^kx_i$, which are collected from all the sources and then optimized by the system coordinator by solving the following problem

$$\begin{aligned} \begin{array}{l} {\mathop {\max }\limits _{x}} \sum \limits _{i \in \mathcal{S}} {{B_i}({x_i})} \\ \sum \limits _{i \in \mathcal{S}(l)} {{x_i} \le 1\quad \forall l \in \mathcal{L}} \\ {x_i} \ge 0 \quad \forall i \in \mathcal{S} \end{array} \end{aligned}$$

(72)

Traditionally utilized first order gradient methods such as the distributed dual gradient and the fast dual gradient method are simulated for comparizon with the proposed DyMonDS-based approach. For the dual gradient algorithm and fast gradient algorithm, the step size for lagrange multiplier increments have been assumed to equal to $\alpha =\frac{2 \sigma }{N_p}{N_s}$ where $\sigma = \nabla ^2 f_i(1)$ is the strong convexity constant for the utility functions in used (Beck et al. 2014). $N_p$ and $N_s$ respectively are the longest path lengths among all sources and maximum number of sources sharing particular link respectively.

For our DyMonDS approach, we do not have to go through the hassle of selecting a step-size however. For all the methods, all of the following termination conditions have been utilized (Beck et al. 2014):

• primal objective function values satisfy $\left| \frac{f(x^{k+1})}{f(x^k)} - 1\right| \le 0.01$

• dual variables satisfy the worst case differences $\Vert \lambda ^{k+1} - \lambda ^{k}\Vert _{\infty } \le 0.01$

• primal feasibility satisfies $A x^k - c < 0.01$ for all the links

Fig. 9

Comparison of the rate of convergence of the proposed DyMonDS approach with that of two other distributed methods for a system with 2 sources and 5 links

All the tested methods converge to a value equal to $x_1 = 0.3$ and $x_2 = 0.7$ for when $C_1 = 10$ and $C_2 = 20$ in the utility functions used in the formulations of NUM in Eq. (69). The convergence rates in terms of the objective function values for all three methods is shown in Fig. 9. Notice the effectiveness of the DyMonDS-based approach which only takes about 10 iterations for convergences in comparison to the other methods which can go upto thousands of iterations. Furthermore, just with a couple of iterations, the approach tends to converge to a near-optimal point which is permissible for large-scale systems.

Similar to the analysis in Beck et al. (2014), we have further produced random networks with a random number of sources in the range [1, 25] and a random number of links in the range of [1, 40] in order to test the scalability. For each of the 50 trials performed, the number of iterations for convergence is shown in the Fig. 10. As anticipated the DyMonDS approach requires very few number of iterations irrespective of the size of the system.

Fig. 10

Comparison of iterations for random network