With the gradual depletion of fossil fuels and escalating environmental issues, traditional energy production and consumption methods are becoming unsustainable. Consequently, large-scale development and utilization of renewable energy have become a strategic focus for sustainable energy development in many countries. China is committed to achieving carbon neutrality by 2060, and high-proportion grid integration of renewable energy will be a key feature of future power systems. As electric cars rapidly develop and find applications, they are increasingly becoming a crucial solution to energy and environmental challenges, emerging as flexible loads in power grids. When connected to the grid, electric cars can exchange power bidirectionally: during charging, they act as consumers, and during discharging, they serve as generators. Supported by national policies, the China EV market has experienced rapid growth, with extensive infrastructure development and incentives such as purchase subsidies to encourage consumer adoption.
The interaction between electric cars and power grids has been widely studied, with some research focusing on vehicle-to-grid (V2G) based frequency regulation. Studies indicate that electric cars possess unique advantages in rapidly responding to frequency fluctuations, enabling grid frequency regulation through reasonable charging and discharging strategies. A major challenge arises from the large frequency fluctuations caused by integrating intermittent renewable energy into isolated grids, as the capacity of load frequency control (LFC) may be insufficient to handle mismatches between generation and demand. LFC is responsible for maintaining power balance in the system, ensuring frequency variations remain within acceptable limits and preserving dynamic performance. In multi-area power systems, LFC is used to balance generation and load demand. However, traditional interconnected power systems often employ a two-layer hierarchical control for economic load dispatch and LFC, neglecting transient economic optimization. Economic model predictive control (EMPC) can achieve economic load scheduling and LFC in a single layer, enabling transient optimization and time-varying control. Yet, due to the arbitrary nature of stage cost functions in EMPC, asymptotic stability cannot be directly proven as in tracking model predictive control (MPC).
To address these issues, we propose a distributed economic model predictive control (DEMPC) method for large-scale electric car charging stations participating in grid auxiliary frequency regulation. First, we establish a model for electric car charging stations assisting in grid LFC. Then, we introduce an EMPC strategy that integrates tracking optimization and economic dispatch into a single layer. Considering multi-area interconnection, we extend EMPC to a distributed control approach, enabling information exchange and collaboration among sub-areas. Finally, we ensure asymptotic stability using an appropriate terminal cost function. Simulation results demonstrate the effectiveness and superiority of this method.
Multi-Area Power Grid Equivalent Control Model
Due to geographical factors, connecting some grid systems to the main grid can be prohibitively expensive. Instead, interconnecting with neighboring grids allows for spatiotemporal resource optimization, significantly reducing operational costs. We consider a four-area interconnected grid, where power exchange between areas is facilitated through AC tie-lines. Each area includes thermal power generation units, consisting of governors, turbines, and power systems. The governor is a fundamental component of thermal power units, and the LFC block diagram for a thermal power unit is shown in the figure, where $u_i$ is the control input from the EMPC controller, $s$ is the complex domain operator, $T_{gi}$ is the governor time constant, $X_{gi}$ is the governor output, $T_{ti}$ is the turbine time constant, $P_{gi}$ is the turbine output power, $P_{di}$ is the load demand, $M_i^a$ and $D_i$ are the equivalent rotational inertia and load damping coefficient, respectively, $\Delta f_i$ is the frequency deviation, and $R_i$ is the speed droop coefficient.
For random load disturbances, the governor acts quickly to achieve primary frequency regulation. The dynamic expressions are as follows:
Governor dynamics:
$$\Delta \dot{X}_{gi} = -\frac{1}{T_{gi} R_i} \Delta f_i – \frac{1}{T_{gi}} \Delta X_{gi} + \frac{1}{T_{gi}} \Delta u_i$$
Turbine dynamics:
$$\Delta \dot{P}_{gi} = -\frac{1}{T_{ti}} \Delta P_{gi} + \frac{1}{T_{ti}} \Delta X_{gi}$$
Generator dynamics:
$$\Delta \dot{f}_i = -\frac{D_i}{M_i^a} \Delta f_i + \frac{1}{M_i^a} \Delta P_{gi} – \frac{1}{M_i^a} \Delta P_{di} – \frac{1}{M_i^a} \Delta P_{\text{tie},i}$$
The tie-line power deviation injected into each area is expressed as:
$$\Delta P_{\text{tie},ij} = \frac{2\pi}{s} \sum_{\substack{j=1 \\ j \neq i}}^M T_{ij} (\Delta f_i – \Delta f_j)$$
where $T_{ij}$ is the synchronization coefficient of the transmission line, and $M$ is the number of sub-areas.
Electric Car Aggregated Charging Station Model
The grid consists of thermal power units, multiple electric cars, and loads (including wind disturbances). Since each electric car charging station contains a different number of electric cars, an equivalent EV model can be considered, parameterized by different inverter capacities for each electric car. The equivalent electric car model for LFC simulates the behavior of an electric car battery, used to compute the total charging/discharging power in a controllable state. The time constant for the electric car is $T_e$, $\Delta u_e$ is the LFC signal sent to the electric car, $\pm \lambda_e$ are the power ramp rate constraints, $\pm \Delta \mu_e$ are the charging and discharging limits, and $\pm \Delta \delta_e$ are the charging/discharging rate limits. $E_{\text{max}}$ and $E_{\text{min}}$ are the maximum and minimum controllable energy of the battery, respectively. $\Delta P_e$ is the charging/discharging power; when $\Delta P_e = 0$, the electric car is in idle mode; when $\Delta P_e > 0$, it is discharging; and when $\Delta P_e < 0$, it is charging. Note that if the energy storage of an electric car exceeds $E_{\text{max}}$, it can no longer charge, and its discharging range is $0$ to $\mu_e$; similarly, if the energy storage falls below $E_{\text{min}}$, the charging range is $-\mu_e$ to $0$. The dynamic model is expressed as:
$$\Delta \dot{P}_{Ek} = \frac{1}{T_e} \Delta u_{Ek} – \frac{1}{T_e} \Delta P_{Ek}$$
Constraints for electric cars in the system:
$$
\begin{aligned}
-\mu_{ek} &\leq \Delta P_{Ek}(t) \leq \mu_{ek}, \quad k = 1, 2 \\
-\delta_{ek} &\leq \Delta \dot{P}_{Ek}(t) \leq \delta_{ek}, \quad k = 1, 2
\end{aligned}
$$
The state of charge (SOC) of an electric car reflects the current operating state and regulation capability of the battery, expressed as:
$$
\begin{aligned}
\text{SOC}_k &= \text{SOC}_0 – \frac{\eta \int P_{ei} dt}{C_n} \\
\Delta \text{SOC}_i &= \text{SOC}_0 – \text{SOC}_k
\end{aligned}
$$
where $\eta$ is the power loss coefficient, $C_n$ is the rated capacity, $\text{SOC}_0$ is the initial state of charge, and $\text{SOC}_k$ is the SOC value at time $k$.
For a 100 MW grid, assuming a 5% frequency regulation capacity is required, each electric car has a rated capacity of 50 kW, an adjustable energy range of 70% to 90%, and a power loss coefficient of 90%. If frequency regulation is entirely handled by electric cars, approximately 560 electric cars would need to be connected to the grid. However, with thermal power units, fewer electric cars are needed to assist in frequency regulation based on the capacity of traditional units.
Grid Model Based on LFC Controller
The proposed LFC controller and grid framework include thermal power units, two equivalent electric car aggregated charging stations (EV1 and EV2), and load disturbances (the load disturbance includes both load and wind power fluctuations). The state-space representation of the model is:
$$
\begin{aligned}
\dot{x}_i(t) &= A_{ii} x_i(t) + B_{ii} u_i(t) + F_{ii} d_i(t) + \sum_{i \neq j} A_{ij} x_j(t) \\
y_i &= C_{ii} x_i(t)
\end{aligned}
$$
where the state vector $x_i = [\Delta f_i, \Delta P_{gi}, \Delta X_{gi}, \Delta P_{E1i}, \Delta P_{E2i}, \Delta \text{SOC}_i, \Delta P_{\text{tie},i}]^T$, $u_i = [u_{gi}, u_{E1i}, u_{E2i}]^T$, $d_i = [\Delta P_{di}]$, and $y_i = [\Delta f_i]$. The matrices $A_{ii}$, $B_{ii}$, $F_{ii}$, and $A_{ij}$ are defined based on the system parameters. Discretizing this model yields:
$$
\begin{aligned}
x_i(k+1) &= A_{ii} x_i(k) + B_{ii} u_i(k) + F_{ii} d_i(k) + \sum_{i \neq j} A_{ij} x_j(k) \\
y_i &= C_{ii} x_i(k)
\end{aligned}
$$
Algorithm Description
Economic Model Predictive Control
To better illustrate the EMPC algorithm, we compare it with standard MPC. Both algorithms use a discrete-time model: $x^+ = f(x, u)$, with system state $x \in X$, control input $u \in U$, state transition function $f(\cdot): X \times U \rightarrow X$, and optimal setpoint $(x_s, u_s)$ provided by the upper layer of hierarchical control.
Standard MPC, or tracking MPC, acts as a control strategy in the lower layer of hierarchical control. The controller aims to minimize the tracking error between the setpoint $(x_s, u_s)$ and the current system state over a finite horizon $N$:
$$
\begin{aligned}
\min_u & \sum_{i=0}^{N-1} \frac{1}{2} \left( \| x(i) – x_s \|_Q^2 + \| u(i) – u_s \|_R^2 \right) + V_f(x(N)) \\
\text{s.t.} & \quad x^+ = f(x, u) \\
& \quad x(i) \in X, u(i) \in U, \quad i \in I_0:(N-1) \\
& \quad x(0) = x, \quad x(N) \in X_f
\end{aligned}
$$
where $N$ is the optimization horizon, $Q$ and $R$ are constant coefficients, $x$ is the initial state, $X_f$ is the terminal region containing $x_s$, and $V_f(\cdot)$ is the terminal penalty. The upper layer economic optimization typically has a longer time scale, while the lower layer setpoint tracking has a shorter time scale.
In EMPC, the economic cost function $l(x, u)$ is directly optimized by the EMPC controller. The EMPC optimization problem is:
$$
\begin{aligned}
\min_{x,u} & \sum_{i=0}^{N-1} l(x(i), u(i)) + V_f(x(N)) \\
\text{s.t.} & \quad x^+ = f(x, u) \\
& \quad x(i) \in X, u(i) \in U, \quad i \in I_0:(N-1) \\
& \quad x(0) = x, \quad x(N) \in X_f
\end{aligned}
$$
Apart from the cost function, the optimization problem is almost identical to standard MPC. Since the economic stage cost function $l(x, u)$ is directly optimized, the proposed EMPC method improves economic performance. EMPC control performance is equivalent to or better than when the system is in steady state.
DEMPC for Interconnected Grids
Considering multi-area interconnection, we propose a DEMPC method. First, we define the control and optimization objectives for the grid, then use DEMPC to solve the optimization problem.
Problem Description: The LFC for interconnected grids is a nonlinear constrained problem. EMPC optimizes rather than tracks an economic cost function. The economic cost function consists of the following components:
- Generation Cost:
$$F_{pi}(k) = \frac{1}{2} a_i \Delta P_{gi}^2(k) + b_i \Delta P_{gi}(k) + c_i$$
where $a_i, b_i, c_i$ are the generation cost coefficients for the thermal plant in area $i$, with $a_i \geq 0$. - Load Frequency Control:
Frequency deviation arises from imbalance between load demand and power supply. When power meets load demand, frequency deviation approaches zero:
$$\Delta f_i \equiv 0 \iff P_{gi} \equiv P_{di} + P_{\text{tie},i} – P_{EV1,i} – P_{EV2,i}$$
The controller penalizes frequency deviations due to supply-demand imbalance to maintain system frequency stability. Thus, the cost function is defined as:
$$F_{ci}(k) = \Delta f_i(k) Q_{ci} \Delta f_i(k)$$
where $Q_{ci}$ is a weighting coefficient, $Q_{ci} \geq 0$. A larger $Q_{ci}$ reduces $\Delta f_i$ faster but may slow convergence to optimal steady state. - Tie-Line Power Control:
The cost function for tie-line power is:
$$F_{di}(k) = \Delta P_{\text{tie},i}(k) Q_{di} \Delta P_{\text{tie},i}(k)$$
where $Q_{di}$ is a weighting coefficient, $Q_{di} \geq 0$. - Electric Car Regulation Cost:
The regulation cost for electric cars is:
$$F_{ei}(k) = a_{ei} \Delta P_{EVi}^2(k) + b_{ei} \Delta P_{EVi}(k)$$
where $a_{ei}$ and $b_{ei}$ are coefficients for battery degradation and electricity price costs, respectively.
Power system constraints are:
$$
\begin{aligned}
P_{gi,\text{min}} &\leq P_{gi} \leq P_{gi,\text{max}} \\
\dot{P}_{gi} &\leq \dot{P}_{gi,\text{max}} \\
\text{SOC}_{i,\text{min}} &\leq \text{SOC}_i \leq \text{SOC}_{i,\text{max}}
\end{aligned}
$$
The first constraint limits generator or turbine power output, the second is the generator rate constraint, and the third is the electric car state of charge constraint.
DEMPC Method: Combining the cost functions, the stage cost function can be written in convex form:
$$l_i(x_i(k), u_i(k)) = F_{pi}(k) + F_{ci}(k) = x_i(k)^T Q_i x_i(k) + x_i(k)^T q_i + c_i$$
where $Q_i \geq 0$. Based on the system model and cost functions, the coefficients are $Q_i = \text{diag}(Q_{ci}, 0.5a_i, 0, a_{ei}, a_{ei}, Q_{di})$ and $q_i = (0, b_i, 0, b_{ei}, b_{ei}, 0)$. Unlike the tracking cost function in standard MPC, the stage cost function in EMPC can take a more general form.
The local cost function for area $i$ considers its own economic optimization objectives, ignoring the impact of control trajectories on other areas. Thus, the local cost function over the prediction horizon is:
$$\phi_i(x, u_i) = \sum_{t=0}^{N-1} l_i(x_i(t|k), u_i(t|k)) + V_{fi}(x_i(N|k))$$
where $u_i(t|k) = [u(0|k)^T, \dots, u(N-1|k)^T]^T$, $x_i(N|k) = [x(0|k)^T, \dots, x(N-1|k)^T]^T$, and $V_{fi}(x_i(N|k))$ is the terminal cost.
To enable collaboration among areas, $\phi_i(x, u_i)$ is replaced by a cooperation-based cost function that considers both local control objectives and coordination with other areas. The distributed cooperation-based EMPC cost function for each sub-area is defined as a strict convex combination of local cost functions:
$$
V(x, u^{*p}_i) = \min_{u_i} \Phi(x, u^{p-1}_1, \dots, u^{p-1}_{i-1}, u_i, u^{p-1}_{i+1}, \dots, u^{p-1}_M) = \min_{u_i} \left( \alpha_i \phi_i(x, u_i) + \sum_{j \neq i} \alpha_j \phi_j(x, u^{p-1}_j) \right)
$$
where $u^{*p}_i$ is the optimal solution, $\alpha_i > 0$, $\sum_{i=1}^M \alpha_i = 1$, and $p$ is the iteration number at time $k$. For the $i$-th subsystem, only the control sequence $u_i$ is computed and updated at iteration $p$ (initialized with $p=0$). $N$ is the number of sub-areas, and each DEMPC optimizes only its sub-area’s control sequence.
The control sequence for sub-area $i$ at discrete time $k$ is defined as:
$$u^p_i(k) = [u^p_i(0|k)^T, u^p_i(1|k)^T, \dots, u^p_i(N-1|k)^T]^T$$
The candidate initial control sequence is:
$$u^0_i(k+1) = [u^p_i(1|k)^T, \dots, u^p_i(N-1|k)^T, K_{N_i} x_{N_i}(N|k)^T]^T$$
where $K_{N_i}$ is the terminal feedback controller.
The optimal solution for sub-area $i$ after $p$ iterations is:
$$u^{*p}_i(k) = [u^{*p}_i(0|k)^T, \dots, u^{*p}_i(N-1|k)^T]^T$$
The optimization problem for sub-area $i$ is:
$$
\begin{aligned}
V(x, u^{*p}_i) &= \min_{u_i} \Phi(x, u^{p-1}_1, \dots, u^{p-1}_{i-1}, u_i, u^{p-1}_{i+1}, \dots, u^{p-1}_M) \\
\text{s.t.} & \quad x_i(k+1) = \sum_{j \in N_i} A_{ij} x_j(k) + B_{ii} u_i(k) + F_{ii} d_i(k) \\
& \quad x(0|k) = x, \quad x(N|k) \in X_f
\end{aligned}
$$
where $X_f$ is the terminal region containing the optimal steady state $x_s$. The requirement $x(N|k) \in X_f$ ensures that $x(N|k)$ must belong to $X_f$. In interconnected grids, power balance must be satisfied at steady state. Thus, power balance is considered in the terminal region constraints. The terminal region condition for each area is:
$$\| P_{gi}(N|k) – P_{di}(N|k) – P_{\text{tie},i}(N|k) + P_{EV1,i}(N|k) + P_{EV2,i}(N|k) \| \leq \epsilon$$
where $\epsilon$ is a constant. If $\epsilon$ is chosen sufficiently small, power balance can be assumed achievable in each area.
The terminal cost function $V_f(\cdot)$ is selected as:
$$V_f(x(N|k)) = x(N|k)^T P_f x(N|k) + x(N|k)^T p_f$$
where $x(N|k) \in X_f$, a sub-level set $\{x | V_f(x) \leq \epsilon_f\}$, and $X_f$ is chosen based on the above condition.
Stability Analysis
EMPC Stability Proof: We make the following assumptions to establish EMPC stability analysis. Assume that the system is strictly dissipative with supply rate $s_i(x_i(t|k), u_i(t|k)) = l_i(x_i(t|k), u_i(t|k)) – l_i(x_{si}, u_{si})$, and there exists a function $\lambda$ satisfying:
$$\lambda_i[x_i((t+1)|k)] – \lambda_i(x_i(t|k)) \leq -\rho_i(x_i(t|k)) + s_i(x_i(t|k), u_i(t|k))$$
where $\rho_i(\cdot)$ is positive.
This assumption is equivalent to the inequality:
$$\min \left( l_i(x_i(t|k), u_i(t|k)) + \lambda_i(x_i(t|k)) – \lambda_i[x_i((t+1)|k)] \right) \geq l_i(x_{si}, u_{si})$$
To establish EMPC stability, the rotated stage cost is defined as:
$$\tilde{l}_i(x_i(t|k), u_i(t|k)) = l_i(x_i(t|k), u_i(t|k)) – l_i(x_{si}, u_{si}) + \lambda_i(x_i(t|k)) – \lambda_i[x_i((t+1)|k)]$$
where $\tilde{l}_i(x_i, u_i) > 0$ for all $x_i, u_i$, and $\tilde{l}_i(x_{si}, u_{si}) = 0$.
The rotated stage cost can be derived as:
$$\tilde{l}(x(t|k), u(t|k)) = \sum_{i=1}^M \alpha_i \tilde{l}_i(x_i(t|k), u_i(t|k))$$
and the terminal cost as:
$$\tilde{V}_f(x(t|k)) = \sum_{i=1}^M \alpha_i \left( x_i(t|k)^T P_{fi} x_i(t|k) + x_i(t|k)^T p_{fi} \right)$$
Assumption 2: There exists a terminal region control law $K_{N_i}$, $u_i = K_{N_i} x_i$, satisfying the inequality:
$$V_{fi}[x_i((t+1)|k)] – V_{fi}(x_i(t|k)) \leq -l_i(x_i(t|k), u_i(t|k)) + l_i(x_{si}, u_{si})$$
Lemma 1: If the above holds, then the following inequality holds:
$$\tilde{V}_{fi}[x_i((t+1)|k)] – \tilde{V}_{fi}(x_i(t|k)) \leq -\tilde{l}_i(x_i(t|k), u_i(t|k)) + \tilde{l}_i(x_{si}, u_{si})$$
where $\tilde{V}_{fi}(x_i(t|k)) = V_{fi}(x_i(t|k)) + \lambda_i(x_i(t|k)) – V_{fi}(x_{si}) – \lambda_i(x_{si})$.
Proof of Lemma 1:
$$
\begin{aligned}
&\tilde{V}_{fi}[x_i((t+1)|k)] – \tilde{V}_{fi}(x_i(t|k)) \\
&= V_{fi}[x_i((t+1)|k)] + \lambda_i[x_i((t+1)|k)] – \lambda_i(x_{si}) – V_{fi}(x_{si}) \\
&\quad – \left( V_{fi}(x_i(t|k)) + \lambda_i(x_i(t|k)) – \lambda_i(x_{si}) – V_{fi}(x_{si}) \right) \\
&\leq -l_i(x_i(t|k), u_i(t|k)) + l_i(x_{si}, u_{si}) + \lambda_i[x_i((t+1)|k)] – \lambda_i(x_{si}) – \lambda_i(x_i(t|k)) + \lambda_i(x_{si}) \\
&= -\tilde{l}_i(x_i(t|k), u_i(t|k)) + \tilde{l}_i(x_{si}, u_{si})
\end{aligned}
$$
Thus, Lemma 1 is proven.
Using the rotated cost, the new auxiliary cost function is:
$$\tilde{\Phi}(x, u) = \sum_{i=1}^M \alpha_i \sum_{t=0}^{N-1} \tilde{l}_i(x_i(t|k), u_i(t|k)) + \sum_{i=1}^M \alpha_i \tilde{V}_{fi}(x_i(N|k))$$
Lemma 2: The original cost $\Phi(x, u)$ and the rotated cost $\tilde{\Phi}(x, u)$ have the same optimal solution.
Proof of Lemma 2:
$$
\begin{aligned}
\tilde{\Phi}(x, u) &= \sum_{i=1}^M \alpha_i \sum_{t=0}^{N-1} \tilde{l}_i(x_i(t|k), u_i(t|k)) + \sum_{i=1}^M \alpha_i \tilde{V}_{fi}(x_i(N|k)) \\
&= \sum_{t=0}^{N-1} \tilde{l}(x(t|k), u(t|k)) + \tilde{V}_f(x(N|k)) \\
&= \sum_{t=0}^{N-1} \left( l(x(t|k), u(t|k)) – l(x_s, u_s) – \lambda[x((t+1)|k)] + \lambda(x(t|k)) \right) + V_f(x(N|k)) + \lambda(x(N|k)) – V_f(x_s) – \lambda(x_s) \\
&= \sum_{t=0}^{N-1} l(x(t|k), u(t|k)) – \sum_{t=0}^{N-1} l(x_s, u_s) + \lambda(x(0|k)) – \lambda(x(N|k)) + V_f(x(N|k)) + \lambda(x(N|k)) – V_f(x_s) – \lambda(x_s) \\
&= \Phi(x, u) – \left( N l(x_s, u_s) + \lambda(x(0|k)) – V_f(x_s) – \lambda(x_s) \right)_{\text{constant}}
\end{aligned}
$$
Thus, the optimal solution for MPC is the same as for EMPC. Therefore, the optimal solution of Equation (35) can be achieved by solving Equation (23).
Theorem 1: If the conditions hold and the cost function and feasible initial control sequence are defined as in Equations (19) and (22), respectively, then the system is asymptotically stable and converges to $x_s$.
Terminal Cost/Set Calculation
- Using Discrete Linear Quadratic Regulator (DLQR): MPC with terminal cost and terminal set can relax terminal constraints by using a terminal set constraint and adding a terminal cost. Here, we select the terminal set and terminal cost such that the terminal cost is a control Lyapunov function with a corresponding terminal controller. The terminal cost $V_f(x) = x^T P_f x + x^T p_f$ uses a linear-quadratic approach. The terminal cost and controller must satisfy:
$$x^T (A + B K)^T P_f (A + B K) x + p_f^T (A + B K) x – x^T P_f x – p_f^T x \leq -x^T Q x – q^T x$$
Since this must hold for all $x \in X_f$, the linear part must be zero. By contradiction, if the $i$-th component of the linear vector is non-zero, set $x_i = \epsilon$, $x_j = 0$ for $j \neq i$, and let $\epsilon \to 0$; the linear part dominates any quadratic term, and the inequality fails. Eliminating the linear part gives the tracking MPC condition for quadratic terms, with an additional equality constraint for linear terms:
$$
\begin{aligned}
(A + B K)^T P_f (A + B K) – P_f &\leq -Q \\
p_f^T (A + B K – I) &= -q^T
\end{aligned}
$$
The DLQR solution $K, P_f$ satisfies the first inequality. Compute the linear part $p_f$ as:
$$p_f^T = -q^T (A + B K – I)^{-1}$$
Since $(A + B K)$ is Hurwitz with small-magnitude eigenvalues, $(A + B K – I)$ is invertible. - Using Linear Matrix Inequality (LMI): Set the first inequality of Equation (38) as an LMI to optimize additional criteria, such as maximizing the terminal set volume. For the quadratic part, the inequality is:
$$(A + B K)^T P (A + B K) – P \leq -Q$$
Define $E = P^{-1}$, multiply both sides by $E$:
$$E (A + B K)^T E^{-1} (A + B K) E – E \leq -E Q E$$
Define $r = K E$ and move all terms to one side:
$$E – (A E + B r)^T E^{-1} (A E + B r) – E Q E \geq 0$$
To obtain an LMI condition in the optimization variables $(E, r)$, rewrite as:
$$E – \begin{bmatrix} A E + B Y \\ Q^{1/2} E \\ R^{1/2} Y \end{bmatrix}^T \begin{bmatrix} E & 0 & 0 \\ 0 & I & 0 \\ 0 & 0 & I \end{bmatrix}^{-1} \begin{bmatrix} A E + B Y \\ Q^{1/2} E \\ R^{1/2} Y \end{bmatrix} \geq 0$$
Using the Schur complement:
$$\begin{bmatrix} E & E A^T + Y^T B^T & E Q^{1/2} & Y^T R^{1/2} \\ A E + B Y & E & 0 & 0 \\ Q^{1/2} E & 0 & I & 0 \\ R^{1/2} Y & 0 & 0 & I \end{bmatrix} \geq 0$$
Minimize $-\log \det E$ and compute a smaller terminal cost based on LMI. - Computing the Terminal Set: Use $X_f = \{x | x^T P_f x \leq \alpha\}$ instead of the level set of the terminal cost $V_f(x)$, as the level set may not be centered at the origin. Note that since the quadratic inequality holds and $Q \geq 0$, this set satisfies the invariance condition $(A + B K) x \in X_f$. Once $P_f$ and $K$ are computed, compute the scalar parameter $\alpha$ to determine the terminal set size such that input and state constraints are satisfied:
$$K x \in U, \quad \forall x \in X_f, \quad X_f \subseteq X$$
Consider polyhedral state and input constraints:
$$X = \{x | H x \leq h\}, \quad U = \{u | L u \leq l\}$$
Maximizing the terminal set under these constraints can be posed as a linear program (LP). - Lemma 3: The terminal set $X_f = \{x | x^T P_f x \leq \alpha_{\text{max}}\}$ satisfying $X_f \subseteq X$ and $K X_f \subseteq U$ can be computed with the following LP:
$$
\begin{aligned}
\max_{\alpha} & \quad \alpha \\
\text{s.t.} & \quad \| P_f^{-1/2} H_i^T \|^2 \alpha \leq h_i^2, \quad i = 1, 2, \dots, n_x \\
& \quad \| P_f^{-1/2} K^T L_i^T \|^2 \alpha \leq l_i^2, \quad i = 1, 2, \dots, n_u
\end{aligned}
$$
Thus, the terminal cost, terminal controller, and terminal set can be computed as an optimization problem.
Simulation Analysis
We simulate and evaluate the proposed method. The four-area interconnected grid has a power capacity of 100 MW per area, with the number of sub-areas $N=15$ and sampling time $T_s=0.1$ s. The synchronization coefficients are $K_{s12}=3.7$, $K_{s13}=4.0$, $K_{s24}=4.6$, $K_{s34}=4.2$, $K_{s14}=3.5$, $K_{s23}=3.2$ p.u. The power source parameters for each area are summarized in Table 1, and the parameters for electric cars participating in auxiliary frequency regulation are in Table 2. We set $\eta=90\%$, $C_n=50$ kW·h, $Q_{ci}=0.4$, $Q_{di}=0.1$, and load demand changes at $t=3$ s.
| Area | $T_{gi}$ | $T_{ti}$ | $M_i^a$ | $D_i$ | $R_i$ | $a_i$ | $b_i$ | $c_i$ | $P_{di}$ | $\Delta P_{di}$ | $V_i$ |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 0.081 | 0.28 | 3.50 | 2.75 | 2.6 | 3.31 | 0.13 | 1.32 | 0.012 | 0.005 | 1.02 |
| 2 | 0.072 | 0.30 | 3.70 | 3.20 | 2.8 | 2.75 | 0.11 | 1.05 | 0.013 | 0.005 | 1.01 |
| 3 | 0.083 | 0.32 | 4.00 | 2.80 | 2.7 | 3.15 | 0.14 | 1.50 | 0.015 | 0.005 | 0.98 |
| 4 | 0.075 | 0.35 | 3.75 | 2.50 | 2.4 | 3.78 | 0.12 | 1.74 | 0.011 | 0.005 | 1.03 |
| EV | $T_{ei}$ | $\delta_{ei}$ | $\mu_{ei}$ | $E_{\text{max}i}$ | $E_{\text{min}i}$ |
|---|---|---|---|---|---|
| EV1 | 1 | 0.02 | 0.01 | 0.95 | 0.80 |
| EV2 | 1 | 0.02 | 0.015 | 0.90 | 0.75 |
At $t=3$ s, the system experiences a load demand change. The frequency deviations for each area under load demand change are shown in the figure. Due to the load demand change, frequency deviations deviate from zero, indicating that generation does not meet load demand initially. However, with increased generation in each area, the method quickly converges to zero.
The output power for each area under load demand change is shown in another figure. The generation $P_{i,j}$ increases to the optimal setpoint, and due to terminal region constraints, the system state rapidly reaches steady state. The control input $P_{gi}$ responds faster than $\Delta f$ because, in the function relationship, $\Delta f$ has an additional first-order inertial环节 compared to $P_{gi}$ for the same input $u_i$. When $P_{gi}$ reaches steady state, $\Delta f$ remains below the equilibrium point. To accelerate convergence to steady state, $P_{gi}$ is adjusted to exceed the steady-state value, resulting in overshoot between 3-4 s. After overshoot, $\Delta f$ gradually reaches steady state at zero. The black solid line represents the optimal generation power for each area under economic load dispatch in two-layer control, showing that the steady state matches the proposed control method.
The tie-line power injected into each area during load demand change is shown in a separate figure. Here, $P_{\text{tie},i} \geq 0$ indicates area $i$ receives power from the tie-line, while $P_{\text{tie},i} < 0$ indicates area $i$ supplies power to the tie-line.
In traditional grid control, when the system suffers from load changes causing frequency deviations, only thermal power units are adjusted to compensate, leading to limitations in unit output, slow regulation efficiency, and low sensitivity. The output power deviation curves for electric cars and the power curves for the thermal power system are shown in additional figures. For area 1, when the system deviates, thermal power units remain the primary devices for compensating frequency fluctuations. However, compared to systems without electric car participation in frequency regulation, idle electric cars can quickly respond to frequency fluctuations, cooperating with thermal power units until system frequency stabilizes at the reference value.
We compare three different control strategies for the system, each with the same optimization range and economic objective function. Load demand changes by 0.005 p.u. at 1 s in each system. The frequency deviation for area 1 under different control methods is shown in a figure. Under all three strategies, frequency deviation gradually converges to zero after the load demand change. Due to the high dimensionality of the centralized model, simulation computation time is significantly slower than distributed and decentralized approaches, making it unsuitable for industrial applications. In the decentralized model and algorithm, inter-area connections are ignored, resulting in the slowest frequency stabilization. Compared to decentralized EMPC, the distributed EMPC algorithm performance is significantly improved.
Conclusion
We proposed a DEMPC method for electric car charging stations participating in auxiliary frequency regulation with thermal power in multi-area interconnected grids. The economic stage cost function is directly formulated and optimized by the EMPC controller. Based on communication and collaboration among subsystems, the DEMPC controller optimizes the overall system objective function to achieve global optimality. This method enables economic load dispatch, LFC control performance, and transient optimization. Using Lyapunov stability analysis theory, we proven the convergence of the DEMPC optimization algorithm. Simulation results also demonstrate the effectiveness of the optimization algorithm. Compared to centralized and decentralized EMPC methods, this approach offers advantages and holds broad prospects for industrial applications, particularly in the context of China EV development and the integration of electric cars into power systems.