-
Clean energy and water are two essential resources for the development of our society. Traditionally, the production of potable water and power generation are thought of as separable problems. However, with the development of technology, they are becoming increasingly correlated. For example, multi-stage flash (MSF) seawater desalination systems are usually coupled with power plants because they can use the wasted energy of used gas (which exits from gas turbine cycles) to produce potable water, while the gas turbine generator needs water to produce energy. With this scheme, it contributes to improving the fuel efficiency of the whole plant [1].
In recent years, a large number of researchers have investigated the energy-water nexus from various aspects. For example, Ref. [2] developed a model to save the energy and water by optimally controlling the lawn irrigation. Ref. [3] implemented an energy optimization model with a dedicated water module to assess an optimal “water-energy” mix. Ref. [4] developed a quantitative, physics-based model of the energy-water nexus to optimize the energy and water systems from an engineering system perspective. Ref. [5] developed a competitive Markov decision process model for the energy water-climate change nexus and the model was solved by a reinforcement learning algorithm. Ref. [6] established a multi-objective optimization model to study the room and realization path of both energy and water conservation under the prerequisite of stable economic development. Ref. [7] developed an integrated model analysis framework and tool to help predict and satisfy water, energy, and food demands based on model inputs representing productions costs, socioeconomic demands, and environmental controls. Although these research efforts are helpful for solving the energy-water nexus issues, most work adopts centralized methods to handle the associated optimization problem. It means that these approaches require a control center to acquire and process all the data from the whole network. As the scale of network becomes larger and larger, the control center can’t bear so much computation burden.
To overcome the weakness of centralized method, various distributed algorithms have been developed for large-scale network optimization problem. For example, Ref. [8] proposed a distributed adaptive droop control method for DC micro-grid to optimize power dispatch. Ref. [9] designed a distributed method based on incremental cost for the conventional economic dispatch problem (EDP). Ref. [10] still adopted the equal incremental cost criterion to achieve the optimal dispatch, but the proposed algorithm can estimate the mismatch between demand and total power in a collective sense. Ref. [11] firstly used the logarithmic barrier function to reformulate the economic dispatch problem with capacity constraints and then employed the consensus approach to design the distributed algorithm. Ref. [12] proposed a distributed algorithm on the basis of projected gradient and finite-time average consensus algorithms for smart grid systems. Ref. [13] used the method of bisection and a consensus-like iterative method to solve the economic dispatch problem (EDP) in a smart grid. Ref. [14] proposed a distributed algorithm based on Laplacian-gradient dynamics for an EDP without and with capacities constraints. Then, Ref. [15] extended the initial strategy and presented the algorithm which could converge to the optimal solution of the dispatch problem starting from any initial power allocation in a distributed manner. Ref. [16] proposed an algorithm which is capable of solving the EDP in a minimum number of time steps instead of asymptotically. Ref. [17] employed the ADMM and finite-time algorithm to design a distributed algorithm for the EDP over the directed graph. Ref. [18] first designed a distributed inexact consensus-based ADMM method for the unconstrained optimization problem. Based on the work[18], Ref. [19] proposed a distributed inexact dual consensus ADMM to solve the EDP with complex structure object function. Ref. [20] provided two classes of continuous-time algorithms to solve the EDP in an initialization free and scalable manner by adopting either projection or differentiated projection method. Ref. [21] developed a distributed algorithm for the EDP with local non-smooth cost function over weight-balanced digraphs. Ref. [22] adopted a projected output feedback method to solve the EDP with uncertain parameters.
In this paper, we present a novel continuous-time distributed algorithm to solve the economic dispatch problem of the energy-water nexus with local inequality constraints in the framework of non-smooth analysis and algebraic graph theory. To the best of our knowledge, this is the first attempt to provide a quantitative solution for the economic dispatch problem of the energy-water nexus. Firstly, based on the exact penalty function method, we transform the original economic dispatch problem into an equivalent problem without local constraints. Then we propose a distributed algorithm to solve the equivalent problem in a distributed manner, which means that each plant or agent of the network only needs to obtain its neighbors’ information. Therefore, the proposed algorithm can be effectively applied to the large-scale water and power network. Compared with many existed algorithms which cannot get the correct optimal solution if the initial condition has an error, the proposed algorithm in this paper allows a group of plants to solve the economic dispatch problem for any initial value. Moreover, in our algorithm, we do not require the agents or plants to share their respective gradient information with their neighbors. In other words, the proposed algorithm can favorably protect plants’ privacy.
Notations:
$R$ and${R^N}$ represent the set of real numbers and real N-dimensional column vector, respectively;${{{1}}_N}$ (or${{{0}}_N}$ ) denotes an N-dimensional column vector whose all elements are${\rm{1}}$ (or$0$ ); for a vector or a matrix${{x}}$ ,${{{x}}^{\rm {T}}}$ represents its transpose, and$\left\| \cdot \right\|$ represents the Euclidean norm of a vector or the corresponding induced norm of a matrix;$\inf (S)$ represents the greatest lower bound of the set S;${{A}} \otimes {{B}}$ denotes the Kronecker product of the matrix${{A}}$ and${{B}}$ . -
In this section, some basic concepts about algebraic graph theory are firstly introduced. Then we review some notions from non-smooth analysis and differential inclusion. Finally, we formulate the economic dispatch problem of the energy water nexus.
-
Now we present some notions from algebraic graph theory[23]. A graph is a triplet
${G = (V;E;{{A}})}$ where V denotes the vertex set,${E \subseteq V \times V}$ represents the edge set, and${{A}}$ is the adjacency matrix. An edge from$j$ to$i$ , denoted by${(i,j)}$ , means that agent$i$ can receive information from agent$j$ . An adjacency matrix is defined by${{{A}} = [{a_{ij}}] \in {R^{N \times N}}}$ , where${{a_{ij}} > {\rm{0}}}$ if${(i,j) \in E}$ and${{a_{ij}} = {\rm{0}}}$ , otherwise. A graph is undirected if${(i,j) \in E}$ and${(j,i) \in E}$ simultaneously. An undirected graph is connected if there is a path between any pair of vertexes. The Laplacian matrix is defined by${{{L}} = {{D}} - {{A}}}$ , where${{ D} = {\rm {diag}}\left( {{d_{\rm{1}}},d_{\rm{2}}\cdot \cdot \cdot,{d_N}} \right)}$ with${d_i} = \displaystyle\sum\limits_{j = 1}^N {{a_{ij}}} $ . For the Laplacian matrix${{L}}$ , at least one of the eigenvalues of${{L}}$ is zero and the rest of them have nonnegative real parts. If$G$ is an undirected and connected graph, then$0$ is a simple eigenvalue of${{L}}$ and the other eigenvalues are positive numbers. Additionally, the eigenvector corresponding to the eigenvalue$0$ is given by${\upsilon {{{1}}_N}}$ for some constant${\upsilon }$ . Throughout this paper, the following assumption is used for the graph${G = (V;E;{{A}})}$ .Assumption 1: The graph
${G = (V;E;{{A}})}$ is undirected and connected. -
We recall some notions from non-smooth analysis[24]. A function
$f:{R^N} \to R$ is locally Lipschitz at${{{x}} \in {R^N}}$ if there exist positive constants$\varepsilon $ and$\delta $ , such that for any vectors${{y}},{{z}} \in B({ x},\delta )$ , one has$|f({{y}}) - f({{z}})| \leqslant \varepsilon \left\| {{{y}} - {{z}}} \right\|$ . If$f$ is Lipschitz near any point${{x}} \in {R^N}$ , then$f$ is said to be locally Lipschitz in${R^N}$ . If the function$f$ is locally Lipschitz in${R^N}$ , then$f$ is differential almost everywhere (a.e.) in the sense of Lebesgue measure. The generalized directional derivative of$f$ at${{x}}$ in the direction${{v}} \in {R^N}$ is defined,$${f^{\rm{0}}}\left( {{{x}};{{v}}} \right) = {\rm {lim}}\;{{\rm {sup}}_{{{y}} \to {{x}},\zeta \to {\rm{0}}}}\frac{{f\left( {{{y}} + \zeta {{v}}} \right) - f({{y}})}}{\zeta }$$ Furthermore, Clarke’s generalized gradient of
$f$ at${{x}}$ is defined by$$\partial f({{x}}) = \{ {{\varsigma}} \in {R^N}:\;{f^0}({{x}},{{v}}) \geqslant ({{\varsigma}} ,{{v}}),\forall {{v}} \in {R^N}\} $$ If
$f:{R^N} \to R$ is a convex function, then one knows that$\partial f({{x}})$ is a nonempty, convex, compact set of${{R^N}}$ , and$\partial f({{x}})$ is upper semicontinuous at${{x}}$ .A time-invariant differential inclusion is given by
$$\dot {{x}}(t) \in F({{x}}(t)),{{x}}(0) = {{{x}}_0},t \geqslant 0 $$ (1) where
$F$ is a set-valued map from${{R^q}}$ to the compact convex subsets of${{R^q}}$ . A solution${{x}}(t):(0,\tau )\to$ ${R^q}$ for${\tau > {\rm{0}}}$ of the differential inclusion (1) is an absolutely continuous curve almost everywhere. A set$M$ is said to be weakly invariant (resp., strongly invariant) with respect to (1), if for every${{{{x}}_{\rm{0}}} \in M}$ ,$M$ contains at least a solution (resp., all the solutions) of (1) starting from${{{x}}_{\rm{0}}}$ . An equilibrium point of (1) is a point${{{{x}}_e} \in {R^q}}$ with${{{{0}}_q} \in F({{{x}}_e})}$ . An equilibrium point${{{z}} \in {R^q}}$ of (1) is Lyapunov stable if, for every${\varepsilon > {\rm{0}}}$ , there exists${\delta = \delta (\varepsilon ) > {\rm{0}}}$ such that, for every initial condition${{{{x}}_{\rm{0}}} = {{x}}({\rm{0}}) \in B({ z};\delta )}$ , every solution${{{x}}(t) \in B({ z};\delta )}$ for all${t > {\rm{0}}}$ .Let
${V{\rm{: }}\,{R^q} \to R}$ be a locally Lipschitz continuous function, and$\partial V$ is the Clarke’s generalized gradient of$V({{x}})$ at${{x}}$ . The set-valued Lie derivative${L_F}V:$ ${R^q} \to R$ of$V$ with respect to the differential inclusion (1) is defined by${L_F}V({ x}) = {\rm{\{ }}a \in R:{{{p}}^{\rm {T}}}{{v}}$ $= a,{{v}} \in F({ x}),$ ${{p}} \in \partial V({ x}){\rm{\} }}$ . In the case when${{L_F}V({{x}})}$ is nonempty, we use${{\rm{max}}{L_F}V({ x})}$ to denote the largest element in the set${{L_F}V({ x})}$ . If$\phi ( \cdot )$ is a solution to (1) and$V:{R^q} \to R$ is locally Lipschitz and regular, then$\dot V(\phi ( \cdot ))$ exists almost everywhere and$\dot V \in {L_F}V(\phi ( \cdot ))$ almost everywhere. In addition, if$V( \cdot )$ is continuous differentiable at${ x}$ , then${L_F}V({{x}}) = \{ {{{v}}^{\rm {T}}}\nabla V({{x}}),v \in F({{x}})\}$ . Next, an invariantce principle is presented for non-smooth regular functions[25].Lemma 1: For the differential inclusion (1), assume that
$F$ is upper semicontinuous and locally bounded, and${F({{x}})}$ takes nonempty, compact, and convex values. Let${V{\rm{: }}\,{R^q} \to\, R}$ be a locally Lipschitz and regular function,${S \subset {R^q}}$ be compact and strongly invariant for (1),$\phi ( \cdot )$ be a solution of (1),$${{X}} = \{ {{x}} \in {R^q}:0 \in {L_F}V({{x}})\} $$ and
$M$ be the largest weakly invariant subset of${{\bar X}} \cap S$ , where${{\bar X}}$ is the closure of X. If${\rm{max}}{L_F}V({{x}}) \leqslant 0$ for all${ x} \in S$ , then$d({{\phi}} (t),M) \to 0$ as$t \to + \infty $ where$d({{\phi}} (t),M) = {\rm{inf}}\{ ||{{\phi}} (t) - {{y}}||,{{y}} \in M\}$ . For scalar$s$ ,${[s]^ + } = s$ if${s > {\rm{0}}}$ , and${[s]^ + } = 0$ , otherwise. -
We now present a mathematical model for co-optimization of power-water nexus. Let
${{x_{{\rm p}_i}} \in R}$ ,${{x_{{\rm w}_j}} \in R}$ denote the power generated by a power plant$i$ and the water produced by a water plant$j$ , respectively. Let${{x_{{\rm cp}_k}} \in R}$ and${{x_{{\rm cw}_k}} \in R}$ denote the power and water produced by a coproduction plant$k$ . Let${{d_{{\rm p}_i}}}$ ,${{d_{{\rm w}_j}}}$ ,${d_{{\rm cp}_k}}$ and${d_{{\rm cw}_k}}$ denote the resource product demands of power plant$i$ , water plant$j$ , and coproduction plant$k$ , respectively. The following notations are introduced to vectorize the formulation:${{{{x}}_{{\rm p}_i}} = {{[{x_{{\rm p}_i}},{\rm{0]}}}^{\rm {T}}}}$ ,${{{{x}}_{{\rm w}_j}} = {{[{\rm{0}},{x_{{\rm w}_j}}]}^{\rm {T}}}}$ ,${{{x}}_{{\rm c}_k}} = {[{x_{{\rm cp}_k}},{x_{{\rm cw}_k}}]^{\rm {T}}}$ ,${{{d}}_{{\rm c}_k}} = {[{d_{{\rm cp}_k}},{d_{{\rm cw}_k}}]^{\rm {T}}}$ ,${{{d}}_{{\rm p}_i}} = {[{d_{{\rm p}_i}},0]^{\rm {T}}},$ ${{{d}}_{{\rm w}_j}} = {[0,{d_{{\rm w}_j}}]^{\rm {T}}}$ . Then the co-optimization problem of power-water nexus is modeled as follows:$$\begin{split} & {{\rm{min}}\sum\limits_{i = {\rm{1}}}^{N_{\rm{p}}} {{f_{{\rm p}_i}}({{{x}}_{{\rm p}_i}}) + \sum\limits_{j = {\rm{1}}}^{N_{\rm{w}}} {{f_{{\rm w}_j}}} } ({{{x}}_{{\rm w}_j}}) + \sum\limits_{k = {\rm{1}}}^{N_{\rm{c}}} {{f_{{\rm c}_k}}({{{x}}_{{\rm c}_k}})} }\\ &\qquad st.\sum\limits_{i = 1}^{{N_{\rm{p}}}} {{{{x}}_{{\rm p}_i}}} + \sum\limits_{j = 1}^{{N_{\rm{w}}}} {{{{x}}_{{\rm w}_j}}} + \sum\limits_{k = 1}^{{N_{\rm{c}}}} {{{{x}}_{{\rm c}_k}}} = \\ &\qquad\;\; \sum\limits_{i = 1}^{{N_{\rm{p}}}} {{{{d}}_{{\rm p}_i}}} + \sum\limits_{j = 1}^{{N_{\rm{w}}}} {{{{d}}_{{\rm w}_j}}} + \sum\limits_{k = 1}^{{N_{\rm{c}}}} {{{{d}}_{{\rm c}_k}}} \\ & {\underline {{x}} _{{\rm p}_i}} \leqslant {{{x}}_{{\rm p}_i}} \leqslant {{{{\bar x}}}_{{\rm p}_i}},\;\;\;\;\;{\rm{for}}\;\;i = 1, 2,\cdot \cdot \cdot,{N_{\rm{p}}}\\ & {\underline {{x}} _{{\rm w}_j}} \leqslant {{{x}}_{{\rm w}_j}} \leqslant {{{{\bar x}}}_{{\rm w}_j}},\;\;\;\;\;{\rm{for}}\;\;j = 1, 2, \cdot \cdot \cdot ,{N_{\rm{w}}}\\ & {\underline {{x}} _{{\rm c}_k}} \leqslant {{{x}}_{{\rm c}_k}} \leqslant {{{{\bar x}}}_{{\rm c}_k}},\;\;\;\;\;{\rm{for}}\;\;k = 1, 2, \cdot \cdot \cdot ,{N_{\rm{c}}} \end{split}$$ (2) where
${f_{{\rm p}_i}}$ ,${f_{{\rm w}_j}}$ and${f_{{\rm c}_k}}$ are respectively the scalar cost functions for the ith power production facility, the jth water production facility and the kth coproduction facility,${{N_{\rm{p}}}}$ ,${{N_{\rm{w}}}}$ and${{N_{\rm{c}}}}$ are the numbers of power, water and coproduction facilities, respectively, and${{{\underline {{x}} }_{{\rm p}_i}}}$ ,${\underline {{x}} _{{\rm w}_j}}$ ,${\underline {{x}} _{{\rm c}_k}}$ ,${\bar {{x}}_{{\rm p}_i}}$ ,${\bar {{x}}_{{\rm w}_j}}$ and${\bar {{x}}_{{\rm c}_k}}$ are positive lower and upper bound vectors, respectively.Remark 1: The equality constraint is a global constraint which denotes the supply-demand balance while the inequality constraints represent the reasonable limits of plants’ production capacity.
Assumption 2: The cost functions
${f_{{\rm p}_i}}$ ,${f_{{\rm w}_j}}$ and${f_{{\rm c}_k}}$ are convex and continuous differentiable. -
Before a distributed algorithm is proposed to solve the optimization problem (2), we need to transform the problem by using the exact penalty function method. We ignore the heterogeneity of the power, water and coproduction plants and define
${{{{{{x}}}}_i} \in {R^{\rm{2}}}}$ as the concatenation vector of${{{x}}_{{\rm{p}}_i}}$ ,${{{x}}_{{\rm{w}}_j}}$ and${{{x}}_{{\rm{c}}_k}}$ , thus the index$i$ ranges from 1 to$N = {N_{\rm{p}}} + {N_{\rm{w}}} + {N_{\rm{c}}}$ . Correspondingly, define${{f_i}({{{x}}_i})}$ as the concatenation vector functions of${{f_{{\rm{p}}_i}}({{{x}}_{{\rm{p}}_i}})}$ ,${{f_{{\rm{w}}_j}}({{{x}}_{{\rm{w}}_j}})}$ ,${f_{{\rm{c}}_k}}({{{x}}_{{\rm{c}}_k}})$ and${{{d}}_i}$ as the demand vector. Based on the optimization problem (2), a transformed model is given as follows:$$ \min {f^\theta }({{x}}) = \sum\limits_{i = 1}^N {f_i^\theta ({{{x}}_i})}\quad {\rm {s.t.}}\;\sum\limits_{i = 1}^N {{{{x}}_i}} = \sum\limits_{i = 1}^N {{{{d}}_i}} $$ (3) where
$$f_i^\theta ({{{x}}_i}) = {f_i}({{{x}}_i}) + \frac{1}{\theta }\sum\limits_{j = 1}^2 {{{(x_i^j - \bar x_i^j)}^ + } + {{(\underline x _i^j - x_i^j)}^ + }} $$ and
$x_i^j$ ,$\underline x _i^j$ and$\bar x_i^j$ denote the jth element of the state variable${{{x}}_i}$ , the lower bound${\underline {{x}} _i}$ and the upper bound${{{\bar x}}_i}$ , respectively. The parameter$\theta $ is a constant and satisfies the following lemma.Lemma 2 [14]: The solutions to the optimization problems (2) and (3) coincide for
${\theta > {\rm{0}}}$ when$$\theta < \frac{1}{{2\mathop {{\rm{max}}}\limits_{x \in {\varDelta _1}} ||\partial f({{x}})|{|_\infty }}}$$ where
$${\varDelta _1}{\rm{ = \{ }}{{x}} \in {R^{2N}}:\sum\limits_{i = 1}^N {{{{x}}_i}} = \sum\limits_{i = 1}^N {{{{d}}_i}} ,{\underline {{x}} _i} \leqslant {{{x}}_i} \leqslant {\bar {{x}}_i},\;\forall \;i\} $$ and the generalized gradient of the function
$f_i^\theta $ with respect to the variable$x_i^j$ is given as follows:$$\frac{{\partial f_i^\theta }}{{\partial x_i^j}}{\rm{ = }}\left\{ {\begin{array}{*{20}{l}} \qquad{\dfrac{{{\rm d}{f_i}}}{{{\rm d}x_i^j}} - \dfrac{{\rm{1}}}{\theta }}&{x_i^j < \underline x _i^j}\\ {\left[\dfrac{{{\rm d}{f_i}}}{{{\rm d}x_i^j}} - \dfrac{{\rm{1}}}{\theta },\dfrac{{{\rm d}{f_i}}}{{{\rm d}x_i^j}}\right]}&{x_i^j = \underline x _i^j}\\ \qquad{\dfrac{{{\rm d}{f_i}}}{{{\rm d}x_i^j}}}&{\underline x _i^j < x_i^j < \bar x_i^j}\\ {\left[\dfrac{{{\rm d}{f_i}}}{{{\rm d}x_i^j}},\dfrac{{d{f_i}}}{{dx_i^j}} + \dfrac{{\rm{1}}}{\theta }\right]}&{x_i^j = \bar x_i^j}\\ \qquad{\dfrac{{{\rm d}{f_i}}}{{{\rm d}x_i^j}} + \dfrac{{\rm{1}}}{\theta }}&{x_i^j > \bar x_i^j} \end{array}} \right.$$ Now, a distributed algorithm is provided to solve the problem (3)
$$\left\{ \begin{array}{l} {{{\dot x}}_i} \in - \partial f_i^\theta ({{{x}}_i}) + {{{\lambda}} _i}\\ {{{\dot \lambda }}_i} = - \displaystyle\sum\limits_{j = {\rm{1}}}^N {({{{\lambda}} _i} - {{{\lambda}} _j}) - \displaystyle\sum\limits_{j = {\rm{1}}}^N {({{{z}}_i} - {{{z}}_j}) + {{{d}}_i} - {{{x}}_i}} } \\ {{{\dot z}}_i} = \displaystyle\sum\limits_{j = {\rm{1}}}^N {({{{\lambda}} _i}} - {{{\lambda}} _j}) \end{array} \right.$$ (4) where
${{{\lambda}} _i} = {(\lambda _i^{\rm{1}},\lambda _i^{\rm{2}})^{\rm {T}}}$ and${{{{z}}_i} = {{({{z}}_i^{\rm{1}},{{z}}_i^{\rm{2}})}^{\rm {T}}}}$ for${i = {\rm{1}}, {\rm{2}},\cdot \cdot \cdot,N}$ . Equivalently, the algorithm (4) can be rewritten as the following vector form,$$\left\{ \begin{array}{l} {{\dot x}} \in - \partial {f^\theta }({{x}}) + {{\lambda}} \\ {{\dot \lambda}} = - ({{L}} \otimes {{{I}}_{\rm{2}}})\lambda - ({{L}} \otimes {{{I}}_{\rm{2}}}){{z}} + {{d}} - {{x}} \\ {{\dot z}} = ({{L}} \otimes {{{I}}_{\rm{2}}}){{\lambda}} \end{array} \right.$$ (5) where
${{{\lambda}} = {{({{\lambda}} _{\rm{1}}^{\rm {T}},{{\lambda}} _{\rm{2}}^{\rm {T}},\cdot \cdot \cdot,{{\lambda}} _N^{\rm {T}})}^{\rm {T}}}}$ ,${{{z}} = {{({{z}}_{\rm{1}}^{\rm {T}},{{z}}_{\rm{2}}^{\rm {T}},\cdot \cdot \cdot,{{z}}_N^{\rm {T}})}^{\rm {T}}}}$ and${{{d}} = {{({{d}}_{\rm{1}}^{\rm {T}}, {{d}}_{\rm{2}}^{\rm {T}}, \cdot \cdot \cdot,{{d}}_N^{\rm {T}})}^{\rm {T}}}}$ .Theorem 1: Under assumption 1 and assumption 2, if
${({{{x}}^*},{{{\lambda}} ^*},{{{z}}^*})}$ is an equilibrium point of the distributed algorithm (5), then${{{{x}}^*}}$ is an optimal solution of the problem (3).Proof: Since
${({{{x}}^*},{{{\lambda}} ^*},{{{z}}^*})}$ is an equilibrium point of the algorithm (5), then one has,$$\left\{ \begin{array}{l} 0 \in - \partial {f^\theta }({{x}}^*) + {\lambda ^*}\\ {\rm{0}} = - ({{L}} \otimes {{{I}}_{\rm{2}}}){\lambda ^*} - ({{L}} \otimes {{{I}}_{\rm{2}}}){{{z}}^*} + {{d}} - {{{x}}^*}\\ {\rm{0}} = ({{L}} \otimes {{{I}}_{\rm{2}}}){{{\lambda}} ^*} \end{array} \right.$$ Because the graph
$G$ is undirected and connected, we have${{{1}}_N^{\rm {T}}{{L}} = {{\textit{0}}}}$ . Furthermore, one has$$\begin{aligned} & \qquad {({{1}}_N^{\rm {T}} \otimes {{{I}}_{\rm{2}}})( - ({{L}} \otimes {{{I}}_{\rm{2}}}){{{\lambda}} ^*} - ({{L}} \otimes {{{I}}_{\rm{2}}}){{{z}}^*} + {{d}} - {{{x}}^*})}=\\ & (({{1}}_N^T{{L}}) \otimes {{{I}}_{\rm{2}}}){{{\lambda}} ^*} - (({{1}}_N^{\rm {T}}{{L}}) \otimes {{{I}}_{\rm{2}}}){{{z}}^*} + ({{1}}_N^{\rm {T}} \otimes {{{I}}_{\rm{2}}})({{d}} - {{{x}}^*})=\\ &\qquad\qquad\qquad\quad {({{1}}_N^{\rm {T}} \otimes {{{I}}_{\rm{2}}})({{d}} - {{{x}}^*}) = {{{0}}_{\rm{2}}}} \end{aligned}$$ Thus one has
$\displaystyle\sum\limits_{i = 1}^N {{{{x}}_i}} = \displaystyle\sum\limits_{i = 1}^N {{{{d}}_i}}$ . Additionally, since${({{L}} \otimes {{{I}}_{\rm{2}}}){{{\lambda}} ^*} = {{{0}}_{{\rm{2}}N}}}$ , one has${{{\lambda}} _{\rm{1}}^* = , {{\lambda}} _{\rm{2}}^* = ,\cdot \cdot \cdot, = {{\lambda}} _N^*}{\rm{ = }}{{\mu}} \in {R^2}$ .Therefore, the equilibrium point
${({{{x}}^*},{{{\lambda}} ^*},{{{z}}^*})}$ of the algorithm (5) satisfies$$ {{0}} \in - \partial {f^\theta }({{{x}}^*}) + {{{1}}_N} \otimes {{\mu}} , \quad \sum\limits_{i = 1}^N {{{{x}}_i}} = \displaystyle\sum\limits_{i = 1}^N {{{{d}}_i}} $$ which is exactly the optimal condition for problem (3). Thus, the equilibrium point
${{{ x}^*}}$ is the optimal solution of the problem (3). The proof is completed.Remark 2: From Theorem 1, it is sufficient to show that the algorithm (5) can converge to the optimal solution
${{{{x}}^*}}$ if the algorithm (5) converges to the equilibrium point${({{{x}}^*},{{{\lambda}} ^*},{{{z}}^*})}$ . Next, we show that the algorithm (5) converges to the equilibrium point.Theorem 2: For any initial state, the solutions of the algorithm (5) converge to the equilibrium point
${({{{x}}^*},{{{\lambda}} ^*},{{{z}}^*})}$ under assumption 1 and assumption 2.Proof: Define an energy function
$$ {V\left( {{{x}},{{\lambda}} ,{{z}}} \right) = \frac{{\rm{1}}}{{\rm{2}}}||{{x}} - {{{x}}^*}|{|^{\rm{2}}} + \frac{{\rm{1}}}{{\rm{2}}}||{{\lambda}} - {{{\lambda}} ^*}|{|^{\rm{2}}}} + \frac{{\rm{1}}}{{\rm{2}}}||{{z}} - {{{z}}^*}|{|^{\rm{2}}} $$ (6) It is noted that there exists a vector
${{{g}} \in \partial {f^\theta }({{x}})}$ such that$$\begin{aligned} & {L_F}V({{x}},{{\lambda}} ,{{z}}) = {({{x}} - {{{x}}^*})^{\rm {T}}}( - {{g}} + {{\lambda}} ) + {({{z}} - {{{z}}^*})^{\rm {T}}}({{L}} \otimes {{{I}}_2}){{\lambda}} +\\ &\quad { {{({{\lambda}} - {{{\lambda}} ^*})}^{\rm {T}}}[ - ({{L}} \otimes {{{I}}_2}){{\lambda}} - ({{L}} \otimes {{{I}}_2}){{z}} + {{d}} - {{x}}]}= \\ & { - {{({{x}} - {{{x}}^*})}^{\rm {T}}}({{g}} - {{{g}}^*}) - {{({{\lambda}} - {{{\lambda}} ^*})}^{\rm {T}}}({{L}} \otimes {{{I}}_2})({{\lambda}} - {{{\lambda}} ^*})} \end{aligned}$$ for any
${{{g}} \in \partial {f^\theta }({{x}})}$ . Since${{f^\theta }({{x}})}$ is convex and the Laplacian matrix is positive semidefinite, one has$$ - {({{x}} - {{{x}}^*})^{\rm {T}}}({{g}} - {{{g}}^*}) - {({{\lambda}} - {{{\lambda}} ^*})^{\rm {T}}}({{L}} \otimes {{{I}}_2})({{\lambda}} - {{{\lambda}} ^*}) \leqslant {\rm{0}}$$ Thus,
${{\rm{max}}{L_F}V({{x}},{{\lambda}} ,{{z}}) \leqslant {\rm{0}}}$ and$V({{x}}(t),{{\lambda}} (t),{{z}}(t))$ is non-increasing. It is easy to obtain that$${V({{x}}(t),{{\lambda}} (t),{{z}}(t)) \leqslant V({{x}}(0),{{\lambda}} (0),{{z}}(0))}$$ Furthermore, the solution
${({{x}}(t),{{\lambda}} (t),{{z}}(t))}$ is bounded. Then we know that$$F({{x}},{{\lambda}} ,{{z}}) = \left( \begin{array}{l} - \partial {f^\theta }({{x}}) + {{\lambda}} \\ - ({{L}} \otimes {{{I}}_2}){{\lambda}} - ({{L}} \otimes {{{I}}_2}){{z}} + {{d}} - {{x}} \\ ({{L}} \otimes {{{I}}_2}){{\lambda}} \end{array} \right)$$ is upper semicontinuous and locally bounded. At the same time, the energy function
${V({{x}},{{\lambda}} ,{{z}})}$ is a locally Lipschitz and regular function. Let$S = {\rm{\{ }}({{x}},{{\lambda}} ,{{z}}):$ $0 \in {L_F}V({{x}},{{\lambda}} ,{{z}}){\rm{\} }}$ . Since${{L_F}V({{x}},{{\lambda}} ,{{z}}) \leqslant {\rm{0}}}$ ,$S$ is a compact and strongly invariant set for the algorithm (5). Thus, according to Lemma 1, the solution of the algorithm (5) starting from any initial value converges to the largest weakly positively invariant set$M$ .Next, we will show that any solution
$({{\bar x}},{{\bar \lambda}} ,{{\bar z}}) \in M$ is an optimal solution of the algorithm (5). Since${{\rm{0}} \in {L_F}V({{\bar x}},{{\bar \lambda}} ,{{\bar z}})}$ , then one has$$\left\{ \begin{array}{l} {({{\lambda}} - {{{\lambda}} ^*})^{\rm {T}}}({{L}} \otimes {{{I}}_2})({{\lambda}} - {{{\lambda}} ^*}) = 0 \\ {({{x}} - {{x}}^*)^{\rm {T}}}({{g}} - {{{g}}^*}) = 0 \end{array} \right.$$ Due to
${\left( {{{L}} \otimes {{{I}}_{\rm{2}}}} \right){{{\lambda}} ^*} = {{{0}}_{{\rm{2}}N}}}$ , one has${\left( {{{L}} \otimes {{{I}}_{\rm{2}}}} \right){{\bar \lambda}} = {{{0}}_{{\rm{2}}N}}}$ . At the same time, it is noted that the Hessian matrix${{\bf H}({{x}})}$ of the function$f({{x}})$ is positive definite, one has$$ {({{x}} - {{{x}}^*})^{\rm {T}}}({{g}} - {{{g}}^*}) = {({{x}} - {{{x}}^*})^{\rm {T}}}{\bf H}(\tau {{x}} + (1 - \tau ){{{x}}^*})({{x}} - {{{x}}^*}) $$ where
${{\rm{0}} \leqslant \tau \leqslant {\rm{1}}}$ . Hence, we have${{{\bar x}} = {{{x}}^*}}$ , which implies that${{\bar x}}$ is an optimal solution of the problem (3). Due to the arbitrariness of$({{\bar x}},{{\bar \lambda}} ,{{\bar z}})$ in$M$ , we can conclude that the equilibrium points in$M$ are all the optimal solutions of the problem (3).Finally, we prove that the solution
${{x}}(t)$ to the algorithm (5) will converge to the largest weakly positively invariant set$M$ . Define${{\phi}} (t) = ({{x}}(t),{{\lambda}} (t),$ ${{z}}(t))$ . Since${\phi (t)}$ is bounded, according to Bolzano-Weierstrass theorem[26], there exists${{{\hat \phi}} = ({{\hat x}},{{\hat \lambda}} ,{{\hat z}})}$ and$({t_k},k = {\rm{1,2}}, \cdot \cdot \cdot )$ such that${{{\phi}} ({t_k}) = ({{x}}({t_k}),{{\lambda}} ({t_k}),}{{z}}({t_k})) \to$ $({{\hat x}},{{\hat \lambda}} ,{{\hat z}})$ as${k \to + \infty }$ . Since${{\rm dist}({{\phi}} (t),M) \to {\rm{0}}}$ , then one has${{{\hat \phi}} \in M}$ and thus${{\hat x}}$ is an optimal solution of the problem (3). Consider a new Lyapunov function$\bar V({{x}}(t),{{\lambda}} (t),$ ${{z}}(t))$ defined as${V({{x}}(t),{{\lambda}} (t),{{z}}(t))}$ by replacing$({{{x}}^*},{{{\lambda}} ^*},{{{z}}^*})$ in the equation (6) with${({{\hat x}},{{\hat \lambda}} ,{{\hat z}})}$ . By similar discussion as above for${V({{x}},{{\lambda}} ,{{z}})}$ , one has${{\rm{max}}{L_F}\bar V({{x}},{{\lambda}} ,{{z}}) \leqslant {\rm{0}}}$ . Due to the continuity of$\bar V$ , for any${\varepsilon > {\rm{0}}}$ , there exists${\omega > {\rm{0}}}$ such that${V({{x}},{{\lambda}} ,{{z}}) \leqslant \varepsilon }$ when$||{{\phi}} - {{\hat \phi}} || \leqslant \omega$ . Since$\bar V$ is monotonically non-increasing on the interval$[0, + \infty )$ , there exists a positive integer$T$ such that when${t \geqslant {t_T}}$ , there is$$\begin{split} & \qquad\qquad\quad {\frac{{\rm{1}}}{{\rm{2}}}||{{x}}(t) - {{\hat x}}|{|^{\rm{2}}}}\leqslant \\ & {\bar V({{x}}(t),{{\lambda}} (t),{{z}}(t)) \leqslant \bar V({{x}}({t_T}),{{\lambda}} ({t_T}),{{z}}({t_T})) < \varepsilon } \\ \end{split} $$ (7) which implies
${{\rm{li}}{{\rm{m}}_{t \to + \infty }}{{x}}\left( t \right) = {{\hat x}}}$ . Furthermore, since$\bar V$ is radially unbounded with respect to${{x}}(t)$ , thus the solution${{x}}(t)$ to the algorithm (5) globally converges to an optimal solution of the problem (3).Remark 3: It is noted that the algorithm (4) is a differential inclusion and implemented by using only the plant’s neighbors information. The gradient information is not shared in the neighborhood and thus the proposed algorithm can favorably protect the privacy of agents.
-
In this section, we apply the algorithm (5) to solve the hybrid optimization problem for an energy-water nexus system. It is assumed that the energy-water nexus system consists of 2 power plants, 2 water plants, and 2 coproduction facilities. The communication topology of the system is an undirected and connected graph which is shown in Fig.1. The cost functions
${f_{{\rm p}_i}}$ ,${f_{{\rm w}_j}}$ ,${f_{{\rm c}_k}}$ of the plants are quadratic functions, i.e.,$$\begin{array}{l} {{f_{{\rm p}_i}} = { x}_{{\rm p}_i}^{\rm {T}}{A_{{\rm p}_i}}{x_{{\rm p}_i}} + {B_{{\rm p}_i}}{{ x}_{{\rm p}_i}} + {C_{{\rm p}_i}}} \\ {{f_{{\rm w}_j}} = { x}_{{\rm w}_j}^{\rm {T}}{A_{{\rm w}_j}}{x_{{\rm w}_j}} + {B_{{\rm w}_j}}{{ x}_{{\rm w}_j}} + {C_{{\rm w}_j}}} \\ {{f_{{\rm c}_k}} = { x}_{{\rm c}_k}^{\rm {T}}{A_{{\rm c}_k}}{x_{{\rm c}_k}} + {B_{{\rm c}_k}}{{ x}_{{\rm c}_k}} + {C_{{\rm c}_k}}} \end{array} $$ The parameters of the plants are given in Table 1. The demand vector is given by
${{{d}} = [{\rm{90,0}},{\rm{20,0}},{\rm{0,90,0}},{\rm{10,50,20,45,40}}}]^{\rm {T}}$ . The parameter$\theta $ in the penalty function${{f^\theta }({{x}})}$ is chosen as${\rm{0}}{\rm{.101\;9}}$ . The initial condition of the algorithm (5) can be randomly given. It is noted that the state${{{{x}}_i} \in {R^{\rm{2}}}}$ of each agents has two components (i.e., power and water), thus we respectively give the evolution of the power and water states of the six plants under the algorithm (5), as shown in Fig. 2 and Fig. 3. In the figures, the black dashed lines denote the optimal solution of the algorithm (5), which is given by${{{x}}_{{\rm p}_{\rm{1}}}^* = {\rm{99}}}$ ,${{{x}}_{{\rm p}_{\rm{2}}}^* = {\rm{80}}}$ ,${{{x}}_{{\rm cp}_{\rm{1}}}^* = {\rm{16}}}$ ,${{{x}}_{{\rm cp}_{\rm{2}}}^* = {\rm{10}}}$ for power and${{{x}}_{{\rm w}_{\rm{1}}}^* = {\rm{50}}}$ ,${{{x}}_{{\rm w}_{\rm{2}}}^* = {\rm{35}}}$ ,${{{x}}_{{\rm cw}_{\rm{1}}}^* = {\rm{40}}}$ ,${{{x}}_{{\rm cw}_{\rm{2}}}^* = {\rm{35}}}$ for water. It can be found that the trajectories of the decision variables about power and water asymptotically converge to the optimal values. Additionally, Fig. 4 and Fig. 5 show that the proposed algorithm can keep the balance of supply and demand of power and water, respectively.Table 1. The parameters of the plants
Plant Type index Max Power Capacity/MW Min Power Capacity/MW Max Water Capacity/ $\rm { {m^3}\cdot h^{-1} }$ Min Water Capacity/ $\rm { {m^3}\cdot h^{-1} }$ power ${i_{\rm{1}}}$ 200 50 0 0 power ${{i_{\rm{2}}}}$ 80 20 0 0 water ${{j_{\rm{1}}}}$ 0 0 50 15 water ${{j_{\rm{2}}}}$ 0 0 35 10 coproduction ${{k_{\rm{1}}}}$ 80 16 40 10 coproduction ${{k_{\rm{2}}}}$ 60 10 75 20 The cost parameters of power plants The cost parameters of water plants ${ {A_{\rm{p}}} }$ ${{B_{\rm{p}}}}$ ${{K_{\rm{p}}}}$ ${{A_{\rm{w}}}}$ ${{B_{\rm{w}}}}$ ${{K_{\rm{w}}}}$ 0.00375 2 5 7 0.00625 1 0.00175 1.75 2 5.8 0.00834 3.25 The cost parameters of coproduction plants ${{A_{{\rm{c}}_{\rm{11}}}}}$ ${{A_{{\rm{c}}_{\rm{12}}}}}$ ${{A_{{\rm{c}}_{\rm{21}}}}}$ ${{A_{{\rm{c}}_{\rm{22}}}}}$ ${{B_{{\rm{c}}_{\rm{1}}}}}$ ${{B_{{\rm{c}}_{\rm{2}}}}}$ ${{K_{\rm{c}}}}$ 0.007433 0.03546 0 0.007093 1.106 4.426 57 0.07881 0.06305 0 0.01261 1.475 5.901 57
Design and Analysis of Distributed Optimization Algorithm for Economic Dispatch Problem of Energy-Water Hybrid Networks
-
摘要: 该文对包含发电设备、产水设备和电水联合生产设备的混合能源网络建立资源经济调度问题的非线性规划模型。针对这类考虑水电生产复杂关系的优化问题,提出一种连续时间分布式算法来寻找经济调度问题的最优解。收敛分析表明,该算法能够在任意初始条件下收敛到最优解,且不需要设备与相邻设备交换成本函数的梯度信息,能够很好地保护设备的隐私信息。最后,数值仿真实验结果验证了求解混合能源网络中资源调度问题算法的性能和有效性。Abstract: In this paper, a mathematical model is firstly built for economic dispatch of a hybrid resource system consisting of power plants, water plants, and coproduction plants. For this kind of optimization model considering the complex relation of the power and water in production, a continuous-time distributed algorithm is proposed to find the optimal solution of the economic dispatch problem. The convergence analysis shows that the trajectories converge to the optimal solution regardless of the initial allocation, avoiding the consideration of the initial procedure in detail. Moreover, applying the proposed algorithm, the plants no longer need to exchange the gradient information of the cost functions with their neighbors, which can protect the privacy of the plants. Finally, a numerical example is provided to demonstrate the performance and effectiveness of the proposed algorithm for the energy-water nexus.
-
Table 1. The parameters of the plants
Plant Type index Max Power Capacity/MW Min Power Capacity/MW Max Water Capacity/ $\rm { {m^3}\cdot h^{-1} }$ Min Water Capacity/ $\rm { {m^3}\cdot h^{-1} }$ power ${i_{\rm{1}}}$ 200 50 0 0 power ${{i_{\rm{2}}}}$ 80 20 0 0 water ${{j_{\rm{1}}}}$ 0 0 50 15 water ${{j_{\rm{2}}}}$ 0 0 35 10 coproduction ${{k_{\rm{1}}}}$ 80 16 40 10 coproduction ${{k_{\rm{2}}}}$ 60 10 75 20 The cost parameters of power plants The cost parameters of water plants ${ {A_{\rm{p}}} }$ ${{B_{\rm{p}}}}$ ${{K_{\rm{p}}}}$ ${{A_{\rm{w}}}}$ ${{B_{\rm{w}}}}$ ${{K_{\rm{w}}}}$ 0.00375 2 5 7 0.00625 1 0.00175 1.75 2 5.8 0.00834 3.25 The cost parameters of coproduction plants ${{A_{{\rm{c}}_{\rm{11}}}}}$ ${{A_{{\rm{c}}_{\rm{12}}}}}$ ${{A_{{\rm{c}}_{\rm{21}}}}}$ ${{A_{{\rm{c}}_{\rm{22}}}}}$ ${{B_{{\rm{c}}_{\rm{1}}}}}$ ${{B_{{\rm{c}}_{\rm{2}}}}}$ ${{K_{\rm{c}}}}$ 0.007433 0.03546 0 0.007093 1.106 4.426 57 0.07881 0.06305 0 0.01261 1.475 5.901 57 -
[1] SHAKIB S E, HOSSEINI S R, AMIDPOUR M, et al. Multiobjective optimization of a cogeneration plant for supplying given amount of power and fresh water[J]. Desalination, 2012, 286: 225-234. doi: 10.1016/j.desal.2011.11.027 [2] WANJIRU E M, XIA X H. Energy-water optimization model incorporating rooftop water harvesting for lawn irrigation[J]. Applied Energy, 2015, 160: 521-531. doi: 10.1016/j.apenergy.2015.09.083 [3] DUBREUIL A, ASSOUMOU E, BOUCKAER T. Water modeling in an energy optimization framework-the water-scarce middle east context[J]. Applied Energy, 2013, 101: 268-279. doi: 10.1016/j.apenergy.2012.06.032 [4] LUBEGA W N, FARID A M. Quantitative engineering systems modeling and analysis of the energy-water nexus[J]. Applied Energy, 2014, 135: 142-157. doi: 10.1016/j.apenergy.2014.07.101 [5] NANDURI V, SAAVEDRA-ANTOLŁNEZ I. A competitive Markov decision process model for the energy-water-climate change nexus[J]. Applied Energy, 2013, 111: 186-198. doi: 10.1016/j.apenergy.2013.04.033 [6] TANG X, JIN Y, FENG C Y, et al. Optimizing the energy and water conservation synergy in china: 2007-2012[J]. Journal of Cleaner Production, 2018, 175: 8-17. doi: 10.1016/j.jclepro.2017.11.100 [7] ZHANG X D, VESSELINOV V V. Integrated modeling approach for optimal management of water, energy and food security nexus[J]. Advances in Water Resources, 2017, 101: 1-10. doi: 10.1016/j.advwatres.2016.12.017 [8] HU J, DUAN J, MA H, et al. Distributed adaptive droop control for optimal power dispatch in dc microgrid[J]. IEEE Transactions on Industrial Electronics, 2018, 65(1): 778-789. doi: 10.1109/TIE.2017.2698425 [9] ZHANG Z A, CHOW M Y. Convergence analysis of the incremental cost consensus algorithm under different communication network topologies in a smart grid[J]. IEEE Transactions on Power Systems, 2012, 27(4): 1761-1768. doi: 10.1109/TPWRS.2012.2188912 [10] YANG S P, TAN S C, XU J X. Consensus based approach for economic dispatch problem in a smart grid[J]. IEEE Transactions on Power Systems, 2013, 28(4): 4416-4426. doi: 10.1109/TPWRS.2013.2271640 [11] LI C J, YU X H, YU W W, et al. Distributed event-triggered scheme for economic dispatch in smart grids[J]. IEEE Transactions on Industrial Informatics, 2016, 12(5): 1775-1785. doi: 10.1109/TII.2015.2479558 [12] GUO F H, WEN C Y, MAO J F, et al. Distributed economic dispatch for smart grids with random wind power[J]. IEEE Transactions on Smart Grid, 2016, 7(3): 1572-1583. doi: 10.1109/TSG.2015.2434831 [13] XING H, MOU Y T, FU M Y, et al. Distributed bisection method for economic power dispatch in smart grid[J]. IEEE Transactions on Power Systems, 2015, 30(6): 3024-3035. doi: 10.1109/TPWRS.2014.2376935 [14] CHERUKURI A, CORTÉS J. Distributed generator coordination for initialization and anytime optimization in economic dispatch[J]. IEEE Transactions on Control of Network Systems, 2015, 2(3): 226-237. doi: 10.1109/TCNS.2015.2399191 [15] CHERUKURI A, CORTÉS J. Initialization-free distributed coordination for economic dispatch under varying loads and generator commitment[J]. Automatica, 2016, 74: 183-193. doi: 10.1016/j.automatica.2016.07.003 [16] YANG T, WU D, SUN Y N, et al. Minimum-time consensus-based approach for power system applications[J]. IEEE Transactions on Industrial Electronics, 2016, 63(2): 1318-1328. doi: 10.1109/TIE.2015.2504050 [17] LI P, HU J P. An ADMM based distributed finite-time algorithm for economic dispatch problems[J]. IEEE Access, 2018, 6: 30969-30976. doi: 10.1109/ACCESS.2018.2837663 [18] JIAN L, ZHAO Y Y, HU J P, et al. Distributed inexact consensus-based ADMM method for multi-agent unconstrained optimization problem[J]. IEEE Access, 2019, 7(1): 79311-79319. [19] JIAN L, HU J P, WANG J, et al. Distributed inexact dual consensus ADMM for network resource allocation[J]. Optimal Control, Applications and Methods, 2019, 40(6): 1071-1087. doi: 10.1002/oca.2538 [20] YI P, HONG Y G, LIU F. Initialization-free distributed algorithms for optimal resource allocation with feasibility constraints and application to economic dispatch of power systems[J]. Automatica, 2016, 74: 259-269. doi: 10.1016/j.automatica.2016.08.007 [21] DENG Z H, LIANG S, HONG Y G. Distributed continuous-time algorithms for resource allocation problems over weight-balanced digraphs[J]. IEEE Transactions on Cybernetics, 2018, 48(11): 3116-3125. doi: 10.1109/TCYB.2017.2759141 [22] ZENG X L, YI P, HONG Y G. Distributed continuous-time algorithm for robust resource allocation problems using output feedback[C]//The 2017 American Control Conference. Seattle, WA, USA: IEEE, 2017: 4643-4648. [23] MESBAHI M, EGERSTEDT M. Graph theoretic methods in multiagent networks[M]. Princeton, New Jersey: Princeton University Press, 2010. [24] AUBIN J P, CELLINA A. Differential inclusions: set-valued maps and viability theory[M]. Berlin: Springer Science & Business Media, 2012. [25] CORTÉS J. Discontinuous dynamical systems[J]. IEEE Control Systems Magazine, 2008, 28(3): 36-73. doi: 10.1109/MCS.2008.919306 [26] BARTLE R G, SHERBERT D R. Introduction to real analysis[M]. Hoboken, NJ: Wiley, 2011.