Nomenclature↩︎

\(C_{\text{buy}}\)

Electricity purchase price.

\(C_{\text{deg}}\)

Battery degradation cost.

\(C_{\text{DG}}\)

Diesel generator cost.

\(C_{\text{sell}}\)

Electricity selling price.

\(E_{\text{max,}i}\)

Battery capacity of the \(i\)th storage.

\(\eta_{\text{batt}}^{\text{ch/dis}}\)

Charging/discharging efficiency of storage.

\(N_{\mathrm{batt}}\)

The number of the battery.

\(N_{\mathrm{br}}\)

The number of the branches.

\(N_{\mathrm{tp}}\)

The number of load types.

\(\Omega_{\text{batt}}\)

Set of battery buses.

\(\Omega_{\text{br/bus}}\)

Set of branches/buses.

\(\Omega_{\mathrm{tp}}\)

Set of load types.

\(P_{\text{buy/sell}}\)

Power purchased from/sold to the main grid.

\(P_{\text{slk}}\)

Power exchanged with the main grid.

\(P_{\text{dis/ch,}i}\)

Discharging/charging power of the \(i\)th storage.

\(P_{\text{batt,}i}\)

Net power exchanged between the \(i\)th energy storage unit and the grid.

\(P_{\text{DG,}i}\)

Power of the \(i\)th diesel generator.

\(P_{\text{load,}i}\)

Predicted load power of the \(i\)th load with.

\(P_{\text{res,}i}\)

Active power of the \(i\)th renewable source.

\(\sigma_{\text{SoC,}i}\)

State of charge (SoC) of the \(i\)th storage.

\(\tilde{l}\)

Auxiliary variable for squared branch current \(|I|^2\).

\(\tilde{v}\)

Auxiliary variable for squared node voltage \(|V|^2\).

1 Introduction↩︎

The large-scale integration of photovoltaic (PV) generation into microgrids enhances system flexibility, reduces carbon emissions, and improves overall economic performance [1], [2]. At the same time, the strong climate sensitivity of PV generation introduces significant uncertainty into system operation. In particular, PV output is highly affected by weather conditions such as cloud coverage, which substantially increases the difficulty of accurate forecasting and calls for energy management systems (EMSs) with strong adaptability to PV power fluctuations. Moreover, the high penetration of renewable energy sources (RESs) with variable output intensifies the operational security challenges of microgrids, especially for networks with non-negligible electrical distances. Under such conditions, the EMS is required to explicitly account for power-flow security constraints to ensure safe and reliable operation.

Addressing the first challenge requires the selection of an appropriate EMS architecture. In general, EMS strategies can be categorized into two complementary paradigms: day-ahead planning and online management. Day-ahead EMS focuses on medium- to long-term scheduling of power generation, load profiles, and electricity prices based on forecasted information, with the primary objective of minimizing overall system operating costs [3]–[5]. Owing to its relatively long decision intervals, day-ahead planning typically imposes modest computational requirements and is well-suited for large-scale power systems with complex network models. As a planning-oriented strategy, it plays a critical role in establishing an economically efficient and operationally feasible baseline, particularly for systems with reasonably predictable operating conditions. However, real-time disturbances and forecast uncertainties may lead to deviations from the planned schedules, especially in microgrids with high penetration of RESs. This motivates the integration of online EMS to complement day-ahead planning by enabling timely corrective actions that enhance both operational reliability and economic performance.

Online EMS approaches can be further categorized into two main types: rule-based methods and optimization-based methods. Rule-based EMSs typically rely on predefined heuristics and operating rules derived from prior experience or engineering knowledge, which makes them simple to implement and computationally efficient. For example, a single-step rule-based strategy was introduced in [6] for online scheduling. A decision tree-based EMS method was proposed in [7]. A rule-based controller incorporating multi-step wind power forecasts was developed in [8] to manage stand-alone microgrids. An adaptive, data-driven, rule-based EMS was proposed in [9] by using a reinforcement learning framework. However, their performance strongly depends on the quality of the predefined rules, and they generally lack the flexibility to handle complex operating conditions, multiple conflicting objectives, and rapidly changing system states. These limitations motivate the adoption of optimization-based EMS frameworks.

Model predictive control (MPC) is the most suitable optimization-based EMS approach for this study due to its receding-horizon structure and its ability to explicitly handle multiple objectives and operational constraints [10]–[12]. Numerous studies have demonstrated the effectiveness of MPC-based EMSs in coordinating distributed energy resources under dynamic operating conditions. For instance, the work in [12] comprehensively reviewed MPC applications in photovoltaic-storage-load scheduling, realizing cost reduction and battery life extension. An MPC-based EMS for vehicle-to-grid charging stations was proposed in [13], enabling seamless transition between grid-connected and islanded modes and reducing the cost. More importantly, online forecast updates within MPC enable adaptive responses to unusual weather conditions, thereby mitigating renewable variability and reducing operating costs [14]. However, the performance of MPC critically depends on the accuracy of power forecasting models, which motivates the development of suitable time-series learning methods to support reliable and adaptive online decision-making.

A wide range of time-series learning methods has been developed for renewable-power and load forecasting. From the forecasting-structure perspective, these methods can be broadly categorized into sequence-to-sequence and autoregressive sequence-to-point paradigms. Sequence-to-sequence methods directly map a historical input sequence to a future prediction sequence and are widely adopted in high-capacity neural-network models, such as CNN-, LSTM-, GRU-, and other attention-based architectures [15]–[18]. These models can capture complex temporal dependencies and implicit temporal patterns when sufficient training data are available. However, their modal adaptation is usually learned inside a single high-capacity end-to-end model, and the learned modes are not explicitly organized as representative PV operating profiles. For the online EMS problem considered in this paper, the available profile-level PV training data are limited, and the test profiles may deviate significantly from the representative training distribution due to abrupt weather changes. Under such data-scarce and distribution-shifted conditions, large black-box sequence models may suffer from unstable extrapolation and overfitting, especially when the online controller must repeatedly update forecasts from short measured input windows.

In contrast, autoregressive sequence-to-point methods perform rolling one-step-ahead prediction using shorter input sequences, making them more suitable for cold-start prediction, limited historical records, and frequent online updates in MPC. Representative approaches include neural-network-based models [19]–[21] and kernel-based models [22], [23]. Neural-network-based methods provide flexible nonlinear approximation but usually require more data to maintain robust generalization. Kernel-based methods, particularly KRR-type models, are more data-efficient and stable in small-sample settings, but conventional isotropic or single-distribution kernel models have limited ability to distinguish heterogeneous feature sensitivities and adapt to multiple PV operating modes. Therefore, a forecasting model that combines sample-efficient kernel learning with explicit modal adaptation is needed for online MPC-based EMS under PV forecast deviations.

More critically, although existing forecasting methods can achieve satisfactory performance when the training and testing profiles follow similar distributions, they remain inadequate for PV trajectories characterized by pronounced randomness and multiple operating modes. A single forecasting model trained around a representative profile may exhibit limited extrapolation capability when confronted with atypical or low-generation PV profiles. This motivates the proposed MMDG-AKRR model, which uses an anisotropic kernel to distinguish the power-sequence and time-stamp features, and further organizes multiple AKRR base models into a modal dictionary. Different from attention-based sequence models that learn latent temporal patterns within one global model, the proposed MMDG-AKRR explicitly embeds several base models trained for representative PV operating profiles into a dictionary. During online MPC operation, the current measured input sequence is compared with the representative modal profiles, and similarity-driven weights are used to combine these local base models. In this way, the proposed method serves as a data-efficient and interpretable online correction mechanism for multi-modal PV trajectories, rather than a large end-to-end deep sequence predictor.

The second challenge, mentioned in the first paragraph, lies in developing an MPC-based EMS framework that explicitly accounts for power-flow security constraints while maintaining a favorable balance between computational efficiency and solution quality. For microgrids with short electrical distances, network losses and voltage drops are typically negligible. In such cases, neglecting the nonlinear power-flow model preserves the linear structure of the MPC optimization problem, thereby ensuring high computational efficiency and reliable solution quality. As a result, many existing MPC-based EMS approaches simplify the network power-flow constraints to a node-injection power balance relationship [24]–[30].

However, for medium- to large-scale microgrids with non-negligible electrical distances, such as campus microgrids, the high penetration of PVs necessitates the explicit consideration of network security constraints to ensure safe operation. Directly embedding nonlinear power-flow equations into the MPC formulation results in a nonconvex optimization problem, which significantly degrades computational tractability and solution quality. To address this issue, a second-order cone programming (SOCP) relaxation technique was proposed in [31], providing a convex formulation that enables the incorporation of power-flow constraints while preserving computational efficiency. Nevertheless, this technique has primarily been applied to day-ahead EMS formulations [32] that are insensitive to real-time computational requirements, and its integration into online MPC-based EMS frameworks remains largely unexplored.

To address the above challenges, this paper proposes a power-flow-constrained online MPC-based EMS framework integrated with the proposed MMDG-AKRR forecasting method. The main contributions are summarized as follows:

• A multi-modal dictionary-guided AKRR forecasting method is developed for online PV trajectory correction. The proposed MMDG-AKRR combines an anisotropic kernel for heterogeneous power-sequence and time-stamp features, a dictionary of AKRR base models trained under representative PV operating modes, and an online similarity-driven weighting strategy based on the measured initial PV sequence. Different from large end-to-end neural sequence predictors, the proposed model is designed as a data-efficient and interpretable correction mechanism for PV profiles with limited training samples and significant weather-driven distribution shifts.

• A forecasting-aware SOCP-constrained online MPC-EMS framework is constructed. The proposed EMS couples online MMDG-AKRR forecast updates with rolling-horizon battery/grid dispatch, while embedding an SOCP-relaxed branch-flow security layer to account for voltage and branch-current constraints. This formulation enables repeated online optimization with improved PV-adaptation capability and tractable computational cost.

• A controlled comparative validation is conducted to clarify the roles of the forecasting and EMS components. The forecasting module is compared with NN, MMDG-NN, MKRR, and AKRR baselines to evaluate the effects of the modal dictionary and anisotropic kernel. The EMS framework is compared with LP-based, SOCP-based, and nonlinear power-flow-constrained day-ahead/MPC schemes on 10-, 18-, and 33-bus microgrids to assess forecasting adaptability, grid-security performance, operating cost, and computational efficiency.

The novelty of this work should therefore be interpreted at two levels. The main methodological novelty lies in the MMDG-AKRR forecasting model, which provides a data-efficient multi-modal PV correction mechanism for online MPC. The EMS-side contribution does not claim MPC or SOCP branch-flow relaxation as standalone new techniques. Instead, it lies in integrating the proposed online forecasting correction with an SOCP-constrained rolling MPC workflow for grid-secure microgrid energy management. This integrated design aims to balance PV-forecasting adaptability, operational security, economic performance, and real-time computational tractability.

The remainder of this paper is organized as follows. Section II formulates the overall problem framework. The forecasting models used for MPC are presented in Section III. Section IV details the relaxed branch flow model and the conic formulation. Section V discusses the case study results, and Section VI concludes the paper.

2 Problem Formulation↩︎

This section first presents a general MPC framework. The proposed MMDG-AKRR-based MPC-EMS is then developed within this framework to provide a clear algorithmic structure.

2.1 General Framework of MPC↩︎

MPC is a receding horizon control strategy. A discrete framework of the MPC problem is given in 1 , with prediction horizon \(N_{\text{pre}}\), prediction step size \(\Delta t_{\text{p}}\), and initial step \(k\in\Omega_{\text{d}}\), where \(\Omega_{\text{d}}=\{0,1,\dots,K_{\text{d}}-1\}\) is the decision step set with decision interval \(\Delta t_{\text{d}}\). \[\tag{1} \begin{align} \min_{\mathbf{u}^{[k\to N_{\text{end}}^{[k]}]}}J&=\sum_{n=k}^{k+N_{\text{pre}}-1}\ell(\mathbf{x}^{[n|k]},\mathbf{y}^{[n|k]},\mathbf{u}^{[n|k]})\\ \text{s.t.: }\:\:\:\:\;&\;\notag\\ \mathbf{x}^{[n+1|k]}&=\mathbf{x}^{[n|k]}+\mathbf{f}_{\text{Dyn}}(\mathbf{x}^{[n|k]},\mathbf{y}^{[n|k]},\mathbf{u}^{[n|k]};\Delta t_{\text{p}})\tag{2}\\ \mathbf{y}^{[n|k]}\;\;\;&=\mathbf{f}_{\text{Alg}}(\mathbf{x}^{[n|k]},\mathbf{y}^{[n|k]},\mathbf{u}^{[n|k]})\tag{3}\\ \;\mathbf{x}^{[n|k]}\in&\Omega_{\mathbf{x}},\;\mathbf{y}^{[n|k]}\in\Omega_{\mathbf{y}},\mathbf{u}^{[n|k]}\in\Omega_{\mathbf{u}},\tag{4} \end{align}\] where the notation \(\Box^{[n|k]}\) denotes the value of a variable at prediction step \(n\), given the initial decision step \(k\). \(N_{\text{end}}^{[k]}=k+N_{\text{pre}}-1\) is the last index of \(\mathbf{u}\). Equation 2 denotes the dynamic constraints \(C_{\text{dyn}}\), including the state of charge (SoC) dynamics and power forecasting models in the EMS setting. Equation 3 represents algebraic constraints \(C_{\text{alg}}\) such as power balance and grid security constraints. Equation 4 collects other inequality constraints \(C_{\text{ieq}}\) on state, output, and control variables.

2.2 Overview of the Proposed MPC-based EMS↩︎

This section presents the detailed optimization model, as shown in Fig. 1, for the proposed MPC-EMS based on the framework in 1 . The decision variables, state variables, and algebraic variables are listed below, respectively, \[\begin{align} \mathbf{u}=&\{{P}_{\text{buy}},{P}_{\text{sell}}, {P}_{\text{dis,}i\in\Omega_{\text{batt}}},{P}_{\text{ch,}i\in\Omega_{\text{batt}}},\,{P}_{\text{DG,}i\in\Omega_{\text{DG}}},\notag\\ &\tilde{v}_{i\in\Omega_{\text{bus}}},\tilde{l}_{ij\in\Omega_{\mathrm{br}}},\,P_{ij\in\Omega_{\text{br}}},\Delta C_{\text{inc,}p\in\{1,\dots,N_{\mathrm{tp}}\}}\},\notag\\ \mathbf{x}=&\{\sigma_{\text{SoC,}i\in\Omega_{\mathrm{batt}}},{P}_{\text{res,}j\in\Omega_{\mathrm{res}}},{P}_{\text{load,}j\in\Omega_{\mathrm{load}}}\},\notag\\ \mathbf{y}=&\{P_{j\in\Omega_{\mathrm{bus}}}\}, \end{align}\] The running cost \(\ell\) is \[\begin{align} \ell&= \Delta t_{\text{d}} \Big( - C_{\text{sell}}^{[n|k]} P_{\text{sell}}^{[n|k]} + C_{\text{buy}}^{[n|k]} \big(P_{\text{buy}}^{[n|k]}+P_{\text{loss}}^{[n|k]}\big) \notag\\ &\quad\;\;+ C_{\text{buy}}^{[n|k]} \!\!\sum_{i\in\Omega_{\text{batt}}}\!\! \big( {P}_{\text{ch,}i}^{[n|k]}(1{-}\eta_{\text{batt}}^{\text{ch}}) + {P}_{\text{dis,}i}^{[n|k]}(1{-}\eta_{\text{batt}}^{\text{dis}}) \big) \notag\\ &\quad\;\;+ C_{\text{deg}}^{[n|k]} \!\!\sum_{i\in\Omega_{\text{batt}}}\!\! \big({P}_{\text{ch,}i}^{[n|k]}+{P}_{\text{dis,}i}^{[n|k]}\big) + C_{\text{DG}}^{[n|k]}\!\sum_{i\in\Omega_{\text{DG}}}\!{P}_{\text{DG,}i}^{[n|k]}\Big) \end{align}\] which mainly considers the main grid selling revenue, purchasing cost, battery degradation cost, charging loss, and diesel generator cost over the prediction horizon. Additionally, the active power loss \(P_{\mathrm{loss}}^{[n|k]}=\sum_{ij\in\Omega_{\mathrm{br}}}\tilde{l}_{ij}^{[n|k]} r_{ij}\) is treated as an extra purchase cost, thereby encouraging branch-loss reduction. Details of the constraints at prediction step \(n\) are listed below.

2.2.1 \(C_{\mathrm{dyn}}\): Dynamic constraints↩︎

These constraints are used to describe the trajectory of SoC, RES power (e.g., PV power), and load power in the prediction horizon. \[\begin{align} &\sigma_{\text{SoC,}i}^{[n+1|k]}=\sigma_{\text{SoC,}i}^{[k|k]}+\sum_{d=k}^{n}\left(\frac{{P}_{\text{ch,}i}^{[d|k]}\,\eta_{\text{batt}}^{\text{ch}}}{E_{\text{max,}i}} - \frac{{P}_{\text{dis,}i}^{[d|k]}}{E_{\text{max,}i}\,\eta_{\text{batt}}^{\text{dis}}}\right)\Delta t_{\text{p}},\tag{5}\\ & {P}_{\text{res,}j}^{[n+1|k]}=f_{\text{R}}({P}_{\text{res,}j}^{[n|k]},P^{[n+1|k]}_{\text{dict,}i^*};\Delta t_{\text{p}}),\tag{6}\\ &{P}_{\text{load,}j}^{[n+1|k]}={P}_{\text{load,}j}^{[n|k]}+f_{\text{L}}({P}_{\text{load,}j}^{[n|k]};\Delta t_{\text{p}}).\tag{7} \end{align}\] The SoC dynamic constraints in 5 is expressed in the summation form to make their linear dependence on decision variables explicit. Equations 6 and 7 represent the forecasting models for renewable sources and loads, respectively. Power forecasting is performed before each MPC optimization iteration and incorporated as a fixed input sequence, thereby preserving the convexity of the resulting MPC problem. Details of the forecasting model are provided in Section 3.

2.2.2 \(C_{\text{alg}}\): Algebraic constraints↩︎

These constraints are composed of the node injection power, power balance, and the power-flow model. \[\begin{align} &{P}_{\text{batt,}j}^{[n|k]}={P}_{\text{dis,}j}^{[n|k]}-{P}_{\text{ch,}j}^{[n|k]},{P}_{\text{dis,}j}^{[n|k]}{P}_{\text{ch,}j}^{[n|k]}=0,\tag{8}\\ & P_{\text{slk}}^{[n|k]}=P_{\text{buy}}^{[n|k]}-P_{\text{sell}}^{[n|k]},P_{\text{buy}}^{[n|k]}P_{\text{sell}}^{[n|k]}=0, \tag{9}\\ & P_{j}^{[n|k]} = f_{\mathrm{grid}}(P_{ij}^{[n|k]},P^{[n|k]}_{jk',k'\in\mathcal{C}_{j}},I_{ij}^{[n|k]},r_{ij}) ,\tag{10}\\ & \tilde{V}_{j}^{[n|k]}=f_{\mathrm{node}}(\tilde{V}_{i}^{[n|k]},I_{ij}^{[n|k]},r_{ij}),\tag{11}\\ & \tilde{I}_{ij}^{[n|k]}=f_{\mathrm{branch}}(P_{ij}^{[n|k]},\tilde{V}_{i}^{[n|k]}).\tag{12} \end{align}\] where 8 represents the battery power constraints and 9 denotes the main grid power exchange constraints. These variable-decoupling strategies facilitate problem convexification. Constraint 10 represents the power balance equation based on the detailed branch flow model, where \(P_{j}=P_{\text{batt,}j}+P_{\text{res,}j}+P_{\text{load,}j}\) denotes the nodal power injection. Constraint 11 corresponds to the nodal voltage equation, and 12 describes the branch current equation. The constraints convexification will be addressed in Section 4. Note that all power variables are defined to be positive when injected into the grid, and their detailed definitions are provided in the Nomenclature section following the abstract.

2.2.3 \(C_{\text{ieq}}\): Inequality constraints↩︎

These constraints consist of power-flow security constraints, state-of-charge (SoC) safety constraints, and SoC state preservation conditions. \[\begin{align} & 0\le\tilde{l}_{ij}^{[n|k]}\le I_{ij,\max}^2,\quad V_{j,\min}^2\le\tilde{v}_{j}^{[n|k]}\le V_{j,\max}^2,\tag{13}\\ & 0\leq P_{\text{buy}}^{[n|k]}\leq P_{\text{sys,max}},~0\leq P_{\text{sell}}^{[n|k]}\leq P_{\text{sys,max}},\tag{14}\\ &0\leq{P}_{\text{dis,}j}^{[n|k]}\leq P_{\text{batt,max}},~0\leq{P}_{\text{ch,}j}^{[n|k]}\leq P_{\text{batt,max}},\tag{15}\\ & P_{\text{DG,min}}\leq{P}_{\text{DG,}j}^{[n|k]}\leq P_{\text{DG,max}},\tag{16}\\ & \sigma_{\text{min}}\leq\sigma_{\text{SoC,}i}^{[n|k]}\leq \sigma_{\text{max}},\tag{17}\\ &\sigma_{\text{SoC,}i}^{[0]}= \sigma_{\text{SoC,}i}^{[K_{\text{d}}]},\tag{18} \end{align}\] where \(I_{ij,\max}\) is the current limit of the branch \(ij\), \(V_{j,\min}\) and \(V_{j,\max}\) are the lower and upper bounds of the \(j\)th bus voltage, \(P_{\text{sys,max}}\) is the power exchange limitation with the main grid, \(P_{\text{batt,max}}\) is the battery power rating, and \(P_{\text{DG,min}}\) combine with \(P_{\text{DG,max}}\) are the boundaries of the DG power. Constraint 17 bounds the SoC in a reasonable range. Constraint 18 promotes day-to-day sustainability by preventing systematic under-charging across days.

3 Forecasting Model↩︎

A kernel-based auto-regressive framework, i.e., kernel ridge regression (KRR), which utilizes the data efficiently and is suitable for a small sample scenario, is introduced first, and then the MMDG-AKRR is developed for solar power forecasting.

3.1 Kernel Ridge Regression↩︎

The auto-regressive forecasting model, initialized at step \(k\) within the MPC algorithm, is \[\begin{align} \label{eq95autoReg} P^{[n+1|k]}=P^{[n|k]}+f(P^{[n-N_{L}+1\to n|k]},t_n), \end{align}\tag{19}\] where \(f\) maps the previous \(N_{L}\)-step power sequence, sampled at interval \(\Delta t_{\text{p}}\), to the next-step increment \(\Delta P\). In this study, \(f\) is implemented using kernel ridge regression (KRR), shown in 20 , as it exhibits robust performance in small-sample settings and provides improved generalization, thereby reducing the risk of overfitting. \[\begin{align} \label{eq95KRRseqForecast} f(\mathbf{x}_n)=\mathbf{k}^\top(\mathbf{x}_n)(\mathbf{K}+\lambda \mathbf{I})^{-1}\mathbf{y}_{\text{s}}, \end{align}\tag{20}\] where \(\mathbf{x}_n=[P^{[n-N_{L}+1\to n|k]},t_n]^\top\in\mathbb{R}^{N_L+1}\), \(t_n\) is the normalized timestamp of the last step for the input power sequence \(P^{[n-N_{L}+1\to n|k]}\), and \(\lambda\) denotes the regularization coefficient. The label vector of samples is \(\mathbf{y}_{\text{s}}=[\Delta P_{1},\dots,\Delta P_{N_{\text{sp}}}]^\top\in\mathbb{R}^{N_{\text{sp}}}\) with sample size \(N_{\text{sp}}\). \(\mathbf{K}\in\mathbb{R}^{N_{\text{sp}}\times N_{\text{sp}}}\) is the Gram matrix associated to the Gaussian kernel entries \[k_{\text{ker}}(\mathbf{x}_i,\mathbf{x}_j)=\exp\!(-\frac{\|x_i - x_j\|_2^2}{2\sigma^2})\] with bandwidth \(\sigma\). As an option, one can further incorporate multiple kernels to construct a linear-combination multi-kernel structure, known as multiple-KRR (MKRR) [22], which enhances model expressiveness theoretically at the cost of increased model complexity. \(\mathbf{k}(\mathbf{x}_n)=[k_{\text{ker}}(\mathbf{x}_1,\mathbf{x}_n),\dots,k_{\text{ker}}(\mathbf{x}_{N_{\text{sp}}},\mathbf{x}_n)]^\top\in\mathbb{R}^{N_{\text{sp}}}\) represents the input kernel vector.

While the above KRR-based methods perform well in capturing regular patterns such as load profiles, renewable generation, particularly solar power, often exhibits highly diverse and multi-modal characteristics (Fig. 2). For highly variable PV generation profiles, forecasting models are typically trained using representative typical-day data that cover the most common operating conditions. As a result, such models often exhibit limited transferability to atypical or abrupt scenarios, which may lead to suboptimal planning decisions in the EMS. Moreover, the time-stamp feature differs from power-sequence features in both units and characteristic scales. Directly combining them as a single input may limit the expressive capability of KRR, because a kernel with a single bandwidth has limited flexibility in handling heterogeneous feature sensitivities. Thus, a more reliable forecasting method is required for the MPC-EMS over the prediction horizon.

3.2 Multi-Modal Dictionary-Guided Anisotropic-KRR↩︎

To address the above issues, this study proposes the MMDG-AKRR forecasting method. The main idea is to improve the representation of heterogeneous PV trajectories by combining multiple modal base models rather than fitting all scenarios with a single averaged forecasting model. Specifically, the nonzero time window of daily PV generation profiles is relatively consistent within the same region, and their overall shape exhibits several stable representative patterns. These characteristics motivate the construction of a base-model dictionary for typical groups of PV profiles, where each base AKRR model represents one typical PV operating modal, and the final prediction is obtained through online similarity-based weighting. Since the prediction capability of a modal-dictionary-guided model depends on the representativeness of the predefined modals, the modal dictionary should cover the dominant PV trajectory patterns expected in the considered application. In general, the modal dictionary can be constructed either by clustering historical daily PV profiles, such as using K-means clustering, or by generating representative modal profiles from a typical-day template. The detailed dictionary-construction procedures are provided in 7. It is noted that this approach intentionally neglects high-frequency details in the predicted profiles. However, such details are not critical for EMS applications that focus on power scheduling at relatively large time scales.

The MMDG-AKRR forecasting model used at each prediction step of the MPC horizon, corresponding to 6 , is formulated as \[\begin{align} P^{[n+1]} = P^{[n]} + \Xi^{[n]} \cdot \boldsymbol{\theta}(\mathbf{P}^{[n]}, t_n). \label{eq95MMDGAKRR} \end{align}\tag{21}\] For clarity, the explicit \(k\)-step MPC notation adopted in 19 is omitted here. The weight vector associated with the \(m\) representative PV profile modals is defined as \(\Xi=[w_1,\dots,w_m]\), and the dictionary output vector is given by \(\boldsymbol{\theta}(\mathbf{x}^{[n]}, t_n)=[f_1,\dots,f_m]^\top\). The input \(\mathbf{P}^{[n]}=[P^{[n-N_{\mathrm{L}}+1\to n]}]\) denotes the initial power sequence of length \(N_{\mathrm{L}}\).

Each component \(f_i\) represents a base KRR model trained for the \(i\)th typical PV modal and equipped with an anisotropic Gaussian kernel defined as \[\begin{align} k_{\mathrm{ker}}\big([\mathbf{P}^{[i]},t_i],[\mathbf{P}^{[j]},t_j]\big) =\exp\!\left( -\frac{\|\mathbf{P}^{[i]}-\mathbf{P}^{[j]}\|_2^2}{2\sigma^2} -\frac{\|t_i-t_j\|_2^2}{2\sigma_t^2} \right), \label{eq95AGK} \end{align}\tag{22}\] where \(\sigma\) and \(\sigma_t\) denote the bandwidth parameters associated with the power-sequence features and the time-stamp feature, respectively.

The AKRR hyperparameters \(\lambda\), \(\sigma\), and \(\sigma_t\) are tuned using the proposed multi-start multi-step (MSMS) validation procedure, detailed in the 8.

The overall procedure of the proposed MMDG-AKRR forecasting algorithm is summarized as follows.

• Step 1: Classify the training power profiles into multiple typical modal groups, and decouple each profile into training sample pairs \(\{[\mathbf{P}^{[i]},t_i],\Delta P_i\}_{i=1}^{N_{\text{sp}}}\).

• Step 2: Train a base anisotropic-KRR (AKRR) model \(f\) for each modal group and assemble all base models into the dictionary \(\boldsymbol{\theta}\).

• Step 3: Prior to each prediction update within the MPC loop, evaluate the similarity between the online measured initial \(N_{\text{L}}\)-step input power sequence and the center profile of each training modal. The similarity between the initial sequence and the \(k\)th modal is defined as \[\begin{align} s_k=\Bigg(\max_{n}\frac{1}{N_{\text{L}}}\sum_{i=0}^{N_{\text{L}}-1} \big(P_x^{[i]}-P_{y_k}^{[n+i]}\big)^2\Bigg)^{-1}, \label{eq95similarity} \end{align}\tag{23}\] where \(n=0,\dots,K_{\text{p}}-N_{\text{L}}+1\), \(K_{\text{p}}\) denotes the prediction horizon length, \(P_x^{[i]}\) is the \(i\)th element of the initial input sequence, and \(P_{y_k}^{[n+i]}\) represents the \((n+i)\)th point of the \(k\)th modal center profile. This metric evaluates the inverse of the minimum mean-square error (MSE) over all sliding segments, thereby identifying the most likely modal associated with the current initial condition. Note, the maximum similarity \(S_{\max}=\max_k s_k\) can be used as an indicator of whether the current input sequence is sufficiently represented by the modal dictionary. A small \(S_{\max}\) indicates that the sequence is far from all stored modal profiles, in which case the prediction should be interpreted as an out-of-dictionary forecast rather than a high-confidence modal match. In practice, a threshold for \(S_{\max}\) can be selected from validation profiles or historical operating data to detect truly novel PV behaviors and trigger dictionary enrichment or online model updating when necessary.

• Step 4: Normalize the similarity scores across all modals to obtain the weighting coefficients \(w_k\) in \(\Xi\) using \[\begin{align} w_k=\frac{e^{\beta(s_k-s_{\max})}}{\sum_{i=1}^m e^{\beta(s_i-s_{\max})}}, \label{eq95normalization} \end{align}\tag{24}\] where \(s_{\max}=\max_i s_i\) and \(\beta\) is a tuning parameter controlling the sharpness of the weight distribution. A larger \(\beta\) enforces stronger concentration on the most similar modal. In this study, \(\beta=0.01\) is adopted to achieve a sufficiently sharp distribution while maintaining numerical stability.

• Step 5: Substitute the dictionary \(\boldsymbol{\theta}\) and the weight vector \(\Xi\) into 21 , and iteratively perform the forecasting update until the end of the prediction horizon.

The overall workflow of the proposed MMDG-AKRR algorithm is illustrated in Fig. 3.

4 Power-Flow-Constrained Convex MPC Formulation↩︎

This section first introduces the simplified SOCP-relaxed power-flow/branch-flow model for the MPC-based EMS. Then, the standard conic-program form is organized.

4.1 SOCP Relaxed Power-Flow Model↩︎

The power-flow model [31], [33], illustrated in Fig. 4, is given below \[\tag{25} \begin{align} P_{j} &= \sum_{k' \in \mathcal{C}_{j}} P_{jk'} - P_{ij} + |\tilde{I}_{ij}| ^2 r_{ij}, \tag{26}\\[1ex] Q_{j} &= \sum_{k' \in \mathcal{C}_{j}} Q_{jk'} - Q_{ij} + |\tilde{I}_{ij}|^2 x_{ij}, \tag{27}\\[1ex] \tilde{V}_j &= \tilde{V}_i - \tilde{I}_{ij} \cdot (r_{ij} + j x_{ij}),\text{ with } \tilde{I}_{ij} =\frac{P_{ij} - j Q_{ij}}{\tilde{V}_i^*},\tag{28} \end{align}\]

The non-convexity is introduced by \(|\tilde{I}_{ij}|^2\) and \(\tilde{I}_{ij}\tilde{V}_{i}^*\) terms. Defining auxiliary variables \(\tilde{v}=|\tilde{V}|^2\) and \(\tilde{l}=|\tilde{I}|^2\), the relation \(P_{ij}^2+Q_{ij}^2= \tilde{v}_i\tilde{l}_{ij}\) holds. To obtain a convex SOCP representation, the above equality is relaxed to \(P_{ij}^2+Q_{ij}^2\leq \tilde{v}_i\tilde{l}_{ij}\), which can be rewritten as a second-order cone constraint \(\|[2P_{ij},2Q_{ij},\tilde{v}_i-\tilde{l}_{ij}]\|_2\leq \tilde{v}_i+\tilde{l}_{ij}\). To verify this, the following metric is proposed to quantify the relaxation gap, \[g_{t_{\mathrm d}} = \max_{(i,j)\in\Omega_{\mathrm{br}}} \frac{\big| P_{ij}(t_{\mathrm d})^2 - v_i(t_{\mathrm d}) \ell_{ij}(t_{\mathrm d}) \big|}{\max\!\left(P_{ij}(t_{\mathrm d})^2,\; v_i(t_{\mathrm d}) \ell_{ij}(t_{\mathrm d})\right)} \times 100\% , \label{eq95relaxGap}\tag{29}\] where \(g_{t_{\mathrm d}}\) denotes the maximum relative relaxation gap among all branches at decision step \(t_{\mathrm d}\), reflecting the worst-case tightness of the SOCP relaxation.

In grid-connected operation, the main grid provides strong voltage and reactive power support, such that the microgrid EMS primarily focuses on active power scheduling for economic optimization. Under this condition, reactive power is typically regulated at lower control layers and is not included as an EMS-level decision variable in this study. Moreover, for distribution-level microgrids with relatively high \(R/X\) ratios, i.e., \(x_{ij} \ll r_{ij}\), the coupling between reactive power and system-wide economic objectives is further reduced in the grid-connected mode. It should be noted that reactive power can still provide additional voltage-regulation capability. Therefore, neglecting reactive-power optimization may lead to a conservative assessment of the available voltage-regulation flexibility rather than an overestimation of voltage-security performance. Thus, to reduce computational complexity while preserving sufficient accuracy for EMS-level active-power scheduling, the relaxed branch power-flow model in 25 is simplified by neglecting reactive power optimization, leading to the following reduced formulation.

If reactive-power dispatch or inverter volt-var control is required, the same MPC-EMS framework can be extended by retaining \(Q_{ij}\) and the full \(P\)-\(Q\) SOCP branch-flow constraints in 25 , at the expense of increased decision dimension and computational burden.

\[\tag{30} \begin{align} P_{j} = \sum_{k' \in \mathcal{C}_{j}} P_{jk'} - P_{ij} + \tilde{l}_{ij} r_{ij}, \tag{31}\\[1ex] \tilde{v}_j = \tilde{v}_i + \tilde{l}_{ij} \cdot r_{ij}^2-2P_{ij}r_{ij},\tag{32}\\ \|[2P_{ij},\tilde{v}_i-\tilde{l}_{ij}]\|_2\leq \tilde{v}_i+\tilde{l}_{ij},\tag{33} \end{align}\] where 31 , 32 , and 33 are used to substitute the constraints 10 , 11 , and 12 in the formulated problem. Meanwhile, it provides the information for the active-power-loss term in the objective function, i.e., \(P_{\mathrm{loss}}^{[n|k]}=\sum_{ij\in\Omega_{\mathrm{br}}}\tilde{l}_{ij}^{[n|k]} r_{ij}\).

4.2 Conic Reformulation of MPC↩︎

Several modeling details must be addressed before organizing the MPC problem into a standard conic form. First, to preserve convexity, the mutual exclusivity of the variable pairs \(P_{\mathrm{dis}}\)/\(P_{\mathrm{ch}}\) and \(P_{\mathrm{buy}}\)/\(P_{\mathrm{sell}}\) in 8 and 9 is not explicitly enforced. This is because the linear objective function naturally drives one variable in each pair to zero at optimality.

Similarly, the terminal SoC equality constraint in 18 can be relaxed to the inequality \(\sigma_{\text{SoC,}i}^{[0]}\leq \sigma_{\text{SoC,}i}^{[K_{\text{d}}]}\). This relaxation prevents systematic under-charging across scheduling horizons, while the equality is typically tight under the positive TOU prices considered in this study because an unnecessarily high terminal SoC would require extra charging cost or reduced selling revenue. For special market conditions where such economic tightness may not hold, such as negative grid-purchase prices, the original terminal equality can be retained without changing the proposed MPC-EMS framework.

After convexification, the resulting constraints are grouped into the following categories:

• Box constraints (\(\mathbf{l}_{\text{bd}},\mathbf{u}_{\text{bd}}\)): 13 , 14 , 15 , 16 ;

• SoC constraints (\(\mathbf{A}_{\text{SoC}},\mathbf{b}_{\text{SoC}}\)): 5 , 17 , 18 ;

• Power-flow constraints (\(\mathbf{A}_{\text{grid}},\mathbf{b}_{\text{grid}}\)): 6 , 7 , 8 , 9 , 31 , 32 ;

• Cone constraints (\(\mathbf{A}_{\text{sc }i},\mathbf{b}_{\text{sc }i},\mathbf{d}_{\text{sc }i},\gamma_{\text{sc }i}\)): 33 .

Aggregating the above constraints yields the following standard SOCP formulation of the MPC problem: \[\begin{align} \label{eq95dayahead95SOCP95stform} &\min_{\mathbf{u}_{\mathrm{st}}}\; J=\mathbf{f}^\top\mathbf{u}_{\mathrm{st}}\\ &\text{s.t.:}\,\notag\\ &\mathbf{l}_{\text{bd}}\leq \mathbf{u}_{\mathrm{st}}\leq \mathbf{u}_{\text{bd}},\notag\\ &\mathbf{A}_{\text{SoC}}\mathbf{u}_{\mathrm{st}}\leq\mathbf{b}_{\text{SoC}},\notag\\ &\mathbf{A}_{\Delta E_{\sum}}\mathbf{u}_{\mathrm{st}}\leq\mathbf{b}_{\Delta E_{\sum}},\notag\\ &\mathbf{A}_{\text{grid}}\mathbf{u}_{\mathrm{st}}=\mathbf{b}_{\text{grid}},\notag\\ &\|\mathbf{A}_{\text{sc }i}\mathbf{u}_{\mathrm{st}}-\mathbf{b}_{\text{sc }i}\|_2\leq \mathbf{d}_{\text{sc }i}^\top \mathbf{u}_{\mathrm{st}}-\gamma_{\text{sc }i}, \, i=1,\dots,N_{\text{pre}} N_{\text{br}}\notag \end{align}\tag{34}\]

4.3 Discussion on SOCP Tightness under Weak-Grid and Islanded Conditions↩︎

The SOCP branch-flow relaxation used in this paper follows the standard branch-flow relaxation framework, where the relaxation can be interpreted in two aspects: the conic relaxation of the squared voltage/current magnitude variables and the recovery of voltage/current phase angles. For the full \(P\)-\(Q\) branch-flow model, the conic relaxation can still provide a convex and efficiently solvable formulation for the magnitude-related variables when the feasible region is nonempty. However, exact recovery of the voltage angles is topology-dependent. In particular, for radial distribution networks, the angle recovery is naturally satisfied under standard branch-flow relaxation assumptions, whereas for meshed networks, the angle condition must be explicitly checked or enforced through additional cycle constraints or convexification techniques.

The present study focuses on a grid-connected EMS setting, where the upstream grid provides the voltage reference and primary voltage support. Under islanded operation or weak-grid connection, this assumption becomes weaker or invalid, and inverter voltage/reactive-power support must be modeled more explicitly. Therefore, applying the proposed formulation to islanded or weak-grid microgrids does not necessarily invalidate the convex SOCP formulation itself, but it requires additional verification of angle recovery, relaxation exactness, and dynamic feasibility. In such cases, the full \(P\)-\(Q\) branch-flow model, inverter voltage/reactive-power constraints, and post-optimization AC power-flow or dynamic validation should be included. This extension is beyond the grid-connected active-power EMS scope of this paper and is left for future work.

5 Numerical Study↩︎

This section first introduces the case setup. The effectiveness of the proposed MMDG-AKRR forecasting method is then evaluated. Finally, the proposed MMDG-AKRR-assisted MPC-EMS framework is validated in terms of operational security, cost reduction, and computational efficiency.

The proposed approach is evaluated on three microgrid systems. Detailed results are presented for the 18-bus system, while the results for the 10- and 33-bus systems are summarized subsequently.

5.1 Case Setup↩︎

5.1.1 Load and Solar Profile↩︎

To evaluate the effectiveness of the proposed MPC-based EMS under renewable generation uncertainty, the system performance is examined using a representative operating day characterized by pronounced solar forecast deviations, which commonly arise under variable irradiance conditions such as cloud coverage. This scenario is selected as an illustrative case to highlight the EMS behavior under challenging but realistic conditions, while similar performance trends are observed across other typical operating days. The normalized solar profile used for the following test is shown in Fig. 5 (a), with data obtained from [34]. The normalized residential load and business load profiles used for the test are from [35], as shown in Fig. 5 (b). In the following cases, these normalized profiles are scaled to their respective rated power values.

5.1.2 Battery Information↩︎

The information on the battery used for the tests is as follows. The initial state of SoC is \(\sigma_{\text{SoC}}^{[0]}=0.3\) with bounds \(\sigma_{\text{max}}=0.9\) and \(\sigma_{\text{min}}=0.2\). The charging and discharging efficiency are both \(\eta_{\text{batt}}^{\text{ch}}=\eta_{\text{batt}}^{\text{dis}}=0.95\). The rated energy duration is \(D_{\text{batt}}=4\)h, and the life cycle is \(N_{\text{cyc}}=5000\). The unit device cost is \(C_{\text{unit}}=300\$\)/kWh. The battery energy capacity is defined as \(E_{\text{max}} = P_{\text{battR}} D_{\text{batt}}\), where \(P_{\text{battR}}\) denotes the rated power. Specifically, \(P_{\text{battR}}\) is set to 0.15 MW for all batteries in the 10- and 18-bus cases, while in the 33-bus case it is set to 0.2 MW for Battery 1 and 0.3 MW for Batteries 2–4. The degradation cost is \(C_{\text{deg}}=0.5C_{\text{unit}}/N_{\text{cyc}}\).

5.1.3 Electricity Price Information↩︎

The electricity price information is provided in Table 1.

Here \(C_{\text{buyRes}}\) and \(C_{\text{buyBus}}\) are the purchase tariffs for residential and business (commercial) loads, \(C_{\text{sell}}\) is the sell price to the main grid, and \(C_{\text{DG}}\) is the diesel generator cost. The time-of-use (TOU) periods are defined as follows: Valley, 0h–8h; Off-peak, 8h–16h and 21h–24h; Peak, 16h–21h.

5.1.4 Simulation Environment↩︎

All simulations are conducted in MATLAB on a desktop computer equipped with an Intel Core i5-12400F CPU. The validation problem in 35 for the kernel-based forecasting models is solved using a GA with a population size of 50 and a maximum of 30 generations, which provides a reasonable trade-off between solution quality and computational cost. The EMS optimization problems are solved using either the linprog or coneprog solver, depending on the specific case configuration. The MPC prediction horizon and the decision update interval are set to \(\Delta t_{\mathrm{p}}=1\) and \(\Delta t_{\mathrm{d}} = 0.5\) time units, which corresponds to one hour in this scenario, with a daily-length prediction horizon. The simulation time step is set to \(\Delta t = \Delta t_{\mathrm{p}}/900\).

5.1.5 Comparison Methods↩︎

The forecasting methods considered for performance comparison are summarized as follows:

• NN: An auto-regressive forecasting framework in 19 equipped with a four-layer (128-64-32-1) neural-network (NN)-based inference model;

• MMDG-NN: The proposed MMDG framework equipped with the above neural-network-based inference model;

• MKRR: The proposed MMDG framework equipped with a conventional MKRR model using Gaussian and linear kernels as mentioned in [36];

• AKRR: The KRR forecasting model in 20 using an anisotropic Gaussian kernel;

• MMDG-AKRR: The proposed forecasting method described in Section 3.2;

The EMS frameworks considered for performance comparison are listed as follows:

• LP day-ahead EMS (LP-DA): A conventional linear-programming (LP) day-ahead EMS based on the AKRR forecasting model, which does not explicitly consider detailed branch-flow modeling and associated network security constraints;

• SOCP day-ahead EMS (SOCP-DA): A power-flow-security-constrained day-ahead EMS based on the AKRR forecasting model, where the detailed branch-flow model is incorporated through an SOCP formulation;

• NP day-ahead EMS (NP-DA): A power-flow-security-constrained day-ahead EMS based on the AKRR forecasting model, which adopts the original nonlinear power-flow model and is solved using an interior-point approach;

• LP MPC-EMS (LP-MPC): An MPC version of the LP-based EMS;

• SOCP MPC-EMS (SOCP-MPC): An MPC-based EMS equipped with the AKRR forecasting model and SOCP-relaxed power-flow constraints;

• Proposed MPC-EMS (Proposed MPC): An MPC-based EMS equipped with the proposed MMDG-AKRR forecasting model and SOCP-relaxed power-flow constraints.

For all EMS evaluations, it is assumed that a reliable one-hour-ahead PV power forecast is available at each decision instant, enabled by real-time meteorological measurements/nowcasting and on-site irradiance sensing.

5.2 Results of Forecasting Methods↩︎

To quantify the forecasting accuracy, the mean absolute error (MAE) and the root mean square error (RMSE) are adopted to evaluate the deviation between the predicted power \(\hat{P}\) and the reference trajectory \(P_{\text{ref}}\), defined as \[\begin{align} e_{\text{MAE}}=\frac{1}{K_{\text{p}}}\sum_{i=0}^{K_{\text{p}}}|\hat{P}^{[i]}-P^{[i]}_{\text{ref}}|,\, e_{\text{RMSE}}=\sqrt{\frac{1}{K_{\text{p}}}\sum_{i=0}^{K_{\text{p}}}(\hat{P}^{[i]}-P^{[i]}_{\text{ref}})^2}.\notag \end{align}\] In this study, multi-start forecasting is conducted by initializing the prediction at different time slots, i.e., \(t_0 = 0, 5\!\sim\!17\) h. For each method, the overall forecasting performance is quantified by the average MAE and RMSE over all initial time instants, denoted as \(\bar{e}_{\text{MAE}}\) and \(\bar{e}_{\text{RMSE}}\), respectively.

To ensure generality, training samples are constructed using an auto-regressive sequence-to-point format based on the typical-day scaling strategy mentioned in 7. Specifically, an input sample consists of the previous (\(N_L\))-step PV power sequence and the normalized time stamp, while the output label is the next-step PV power increment. For the MMDG-based models, a modal dictionary is further constructed from scaled versions of the typical profile. In this study, the modal center profiles are generated as (\(a_k P_{\mathrm{typ}}\)), where (\(a_k\in{1.0,0.9,\ldots,0.1}\)). Around each modal center, three nearby profiles, (\(0.96a_k P_{\mathrm{typ}}\)), (\(a_k P_{\mathrm{typ}}\)), and (\(1.04a_k P_{\mathrm{typ}}\)), are used to form the training set for the corresponding base forecasting model. During online prediction, the initial measured PV sequence is compared with these modal center profiles, and the resulting similarity scores are used to determine the weights of the base models.

5.2.1 Regular Day Forecasting↩︎

In this case, the test PV power profile is assumed to lie within the training data distribution. The multi-start forecasting results initialized at different time instants are illustrated in Fig. 6, and the corresponding averaged error metrics are summarized in Table 2.

As shown in Fig. 6, when the forecasting day is consistent with the training data distribution, all methods achieve relatively low prediction errors. Although MKRR exhibits higher error than the other approaches, it still captures the overall trend of the PV profile, indicating that the forecasting models are properly configured. The inferior performance of MKRR can be attributed to its single-bandwidth kernel, which provides limited flexibility in handling heterogeneous features such as the power sequence and time-stamp information, thereby highlighting the advantage of the anisotropic kernel adopted in the AKRR.

Moreover, based on Table 2, the MMDG-based methods exhibit identical performance to their corresponding baseline models. This is because, under regular operating conditions, the initial input sequence shows similar similarity to all modal profiles, leading the MMDG mechanism to assign a dominant weight to the representative typical-day model. As a result, the proposed MMDG framework remains fully compatible with conventional forecasting tasks under regular conditions.

5.2.2 Irregular Day Forecasting↩︎

In contrast, this case considers an irregular PV power profile that lies outside the training data distribution, as illustrated in Fig. 5 (a). The corresponding multi-start forecasting results are presented in Fig. 7, and the averaged error metrics are summarized in Table 3.

As shown in Fig. 7, when no reliable prior information is available for day-ahead forecasting, all methods initialize the prediction at \(t_0 = 0\) h using models trained on typical modal profiles. Consequently, all approaches exhibit noticeable deviations from the true PV trajectory at the beginning of the day, particularly under abrupt and irregular operating conditions. As the MPC prediction horizon recedes and newly measured data are incorporated into the initial input sequence, the forecasting performance of different methods begins to diverge, reflecting their varying abilities to adapt to real-time changes.

Comparing the NN-based results in Tables 2 and 3, it is evident that NN suffers a severe degradation in forecasting accuracy when the test profile deviates from the training distribution. Although AKRR also exhibits increased errors in the irregular-day scenario, its performance remains substantially better than that of NN, indicating stronger generalization capability and robustness to distribution shifts. Furthermore, both NN and AKRR benefit significantly from integration with the proposed MMDG framework, as reflected by the notable error reductions in Table 3. This demonstrates that MMDG effectively enhances model adaptability by exploiting online updated information.

5.2.3 Statistical Forecasting Evaluation on Multiple Irregular PV Profiles↩︎

The preceding irregular-day case provides a representative visualization of the forecasting behavior under a challenging PV trajectory. To further reduce the dependence on a single test profile, an additional statistical evaluation is conducted on multiple irregular PV profiles selected from the same dataset. These profiles cover different generation levels and shape deviations from the typical regular-day profile.

To characterize the selected irregular PV profiles, three profile-level indicators are introduced. Let \(P_{\mathrm{test}}(t)\) denote a selected irregular PV profile and \(P_{\mathrm{typ}}(t)\) denote the typical regular-day profile used for constructing the modal dictionary. Both profiles are normalized by the peak value of the typical profile. The daily energy ratio is defined as \(EnergyRatio = \frac{\sum_{t=1}^{24} P_{\mathrm{test}}(t)}{\sum_{t=1}^{24} P_{\mathrm{typ}}(t)}\). This indicator reflects the overall generation level of the selected irregular profile relative to the typical regular-day profile. The shape deviation is defined as ShapeDeviation = \(\frac{1}{24}\sum_{t=1}^{24} \left| \frac{P_{\mathrm{test}}(t)}{\max_t P_{\mathrm{test}}(t)+\epsilon} \frac{P_{\mathrm{typ}}(t)}{\max_t P_{\mathrm{typ}}(t)+\epsilon} \right|\), where \(\epsilon\) is a small positive constant used to avoid division by zero. This indicator measures the profile-shape mismatch after removing the dominant amplitude effect. The ramping variation is defined as \(RampMean = \frac{1}{23}\sum_{t=1}^{23} \left|P_{\mathrm{test}}(t+1)-P_{\mathrm{test}}(t)\right|\), which describes the average hourly fluctuation intensity of the PV trajectory. The characteristics of the selected irregular PV test profiles are summarized in Table 4.

The statistical forecasting results are illustrated in Fig. 8. For each method, the bar height represents the mean value of the forecasting error over the selected irregular PV profiles. The lower and upper error bars indicate the deviations from the mean to the best-case and worst-case profile-wise errors, respectively, thereby reflecting the performance fluctuation across different irregular profiles. As shown in Fig. 8, MMDG-AKRR achieves the lowest averaged MAE and RMSE among all compared methods. Compared with AKRR, the proposed MMDG-AKRR further reduces both the mean forecasting error and the error fluctuation, indicating that the modal weighting mechanism improves the adaptability of AKRR under irregular PV trajectories. In addition, MMDG-NN also significantly improves the performance of the original NN baseline, although its accuracy remains comparable to the kernel-based baselines. These results suggest that the proposed modal-dictionary strategy is effective for reducing forecasting degradation when the tested PV trajectories deviate from the regular training profile.

These results further verify that the advantage of the proposed MMDG-AKRR model is not limited to a single irregular PV trajectory. Therefore, the proposed forecasting model is adopted in the subsequent MPC-EMS study to provide improved online PV trajectory correction under forecast deviations.

Overall, MMDG-AKRR consistently achieves the lowest forecasting errors among all considered methods. Combined with the better generalization of AKRR over NN and conventional kernel-based approaches, these results support the adoption of MMDG-AKRR within the MPC-EMS framework.

5.3 Results of EMS for the 18-Bus Microgrid↩︎

This case is conducted on the CIGRE 18-bus low-voltage distribution benchmark microgrid system [37], shown in Fig. 9. The branch parameters are determined in accordance with [38]. The power factor of the diesel generator is 0.85. The rated power for three batteries is 0.15MW. Solar panels connected to bus 15 and bus 16 have a rated power of 0.2MW and 0.3MW, respectively. To mitigate voltage drops caused by resistive losses along the long main feeder, a 0.1MW DG is installed at Bus 10, operating within a \(20\%\text{–}100\%\) rated power range. Aggregated residential users on each bus are assigned rated powers \(P_{\mathrm{res}} = (10+20X)\) kW, where \(X \sim \mathcal{U}(0,1)\). The rated power for the business users on buses 2 and 17 is 60 and 100kW, respectively.

5.3.1 Grid safety results↩︎

The bus voltage profiles for all six EMS cases are illustrated in Fig. 10 (a), where the dashed line indicates the voltage security limit. The day-ahead LP-based EMS exhibits pronounced under-voltage during midday, whereas the LP MPC-based EMS suffers severe under-voltage issues in the morning and late afternoon. In contrast, the EMS schemes incorporating power-flow security constraints substantially reduce both the duration and magnitude of voltage violations.

The branch current profiles of the most heavily loaded branch are shown in Fig. 10 (b), where the dashed line denotes the maximum allowable current magnitude. As observed, both the LP day-ahead and LP MPC-based EMS cases exhibit large and sustained overcurrent violations, whereas methods constrained by power-flow security significantly suppress the current overlimit. In particular, the proposed MPC only shows a very small and short-duration current excursion under the strict pointwise threshold. Such minor post-simulation excursions are mainly caused by inverter tracking dynamics, interpolation, and finite simulation-step effects after the optimized dispatch is applied to the dynamic microgrid model, and they can be further mitigated by using a denser decision/simulation step in deployment. Overall, these results demonstrate the necessity of incorporating detailed grid-security constraints into EMS design.

The detailed quantitative overlimit statistics are further reported in Appendix 9.

5.3.2 Operation results↩︎

The battery SoC trajectories for all considered cases are illustrated in Fig. 11. As shown, the SoC of each battery remains within the prescribed bounds of \([0.2,0.9]\). For the convex EMS cases, the batteries are recharged close to their initial SoC levels at the end of the day, which is consistent with the relaxed terminal SoC requirement that the final SoC should be no lower than the initial SoC. The only exception is the NP-DA benchmark. This abnormal terminal-SoC behavior should not be interpreted as looseness of the terminal-SoC relaxation at the optimum. Instead, NP-DA directly embeds the original nonlinear power-flow constraints into the EMS problem, which destroys convexity and makes the interior-point solver unable to guarantee global optimality under practical solver tolerances. Therefore, the NP-DA result reflects the numerical difficulty of directly solving the nonlinear power-flow-constrained EMS problem and further motivates the proposed convex SOCP-based formulation. Compared with the two LP-based cases, the remaining methods charge the batteries more gradually over an extended time interval, thereby mitigating potential grid security issues. Moreover, for the three MPC-based EMS cases, the charging strategy is adaptively rescheduled in response to solar power shortages, avoiding unnecessary full charging of all batteries during non-valley periods and thus reducing additional operational costs.

Aggregated power profiles for each microgrid component are illustrated in Fig. 12. All three day-ahead EMS schemes, namely LP-DA, SOCP-DA, and NP-DA, closely follow the day-ahead scheduling results. This leads to pronounced battery charging at noon (orange curves), which results in suboptimal operating costs, as discussed previously. For the LP-DA case, where power-flow constraints are not considered, such aggressive midday charging further induces branch overcurrent and bus under-voltage issues, as evidenced in Fig. 10 (a) and Fig. 10 (b).

In contrast, the three MPC-based EMS methods—LP-MPC, SOCP-MPC, and the proposed MPC—can adaptively reschedule the battery charging strategy in response to real-time variations in solar power. In particular, these MPC schemes proactively increase battery charging during the valley-price period before anticipated solar power shortages, thereby reducing overall operating costs. However, without explicit power-flow constraints, the LP-MPC scheme still exhibits overly aggressive pre-noon charging behavior, which again leads to grid security violations. By comparison, the two power-flow-constrained MPC schemes distribute the charging power more smoothly over time, ensuring secure operation.

Furthermore, when comparing SOCP-MPC with the proposed MPC, the latter benefits from the improved forecasting accuracy provided by the MMDG-based prediction model. This enables the proposed MPC to more effectively leverage online information, shifting a greater portion of battery charging from the non-valley noon period to the early morning hours, thereby further reducing operating costs while maintaining grid security.

5.3.3 Operation economic and time cost analysis↩︎

Table 5 summarizes both the operating cost and the average decision-time cost per iteration, accounting for the entire prediction horizon. By comparing the day-ahead EMS and the corresponding MPC-based EMS schemes (i.e., LP-DA vs. LP-MPC and SOCP-DA vs. SOCP-MPC), it can be observed that the MPC-EMS methods consistently achieve lower operating costs, since they can reschedule control actions online using updated system information.

Comparing power-flow-constrained methods with their non-constrained counterparts (i.e., LP-DA vs. SOCP-DA and LP-MPC vs. SOCP-MPC), the power-flow-constrained approaches yield lower costs. This improvement arises because EMS schemes with explicit power-flow models can exploit branch-level loss information and thereby enable more efficient power-flow optimization.

The proposed MPC-EMS method further reduces the operating cost by incorporating a forecasting model with higher sensitivity to modal variations during online updates, which leads to more accurate power predictions and improved scheduling decisions.

From a computational perspective, LP-based methods (LP-DA and LP-MPC) exhibit the highest time efficiency. The convex model-based methods, including SOCP-DA, SOCP-MPC, and the proposed MPC approach, also maintain relatively low computational costs and remain suitable for real-time implementation. In contrast, NP-DA, which relies on a nonlinear power-flow model, incurs a significantly higher computational burden and may suffer from suboptimal solutions. This comparison highlights the necessity and effectiveness of cone-relaxed power-flow constraints in achieving a favorable balance between optimality and computational efficiency.

Overall, the proposed MMDG-AKRR-assisted MPC-EMS method achieves the best overall performance by effectively reducing operating costs while ensuring grid security and maintaining time-efficient operation.

5.3.4 Prediction-Horizon Sensitivity Analysis↩︎

The performance of MPC-based EMS may depend on the selected prediction horizon because the controller optimizes the current action based on the predicted future PV generation, load demand, and electricity price. To examine the influence of this factor, an additional prediction-horizon sensitivity test is conducted. Three representative horizons are considered, namely 6 h, 12 h, and 24 h, corresponding to short-, medium-, and long-horizon settings, respectively. The decision update interval remains unchanged, and the controller is rolled forward over the complete 24 h operation period.

The baseline SOCP-MPC and the proposed MPC are compared under each prediction horizon. The total cost, average decision time, power-flow constraint status, and terminal SoC constraint status are summarized in Table 6. For the 6 h and 12 h horizons, the hard end-of-day terminal SoC constraint is not imposed inside every short receding window to avoid artificial infeasibility. Instead, the terminal SoC deficit after the complete 24 h operation is converted into an equivalent charging compensation penalty and included in the reported total cost. Specifically, the missing battery energy is compensated using the peak grid charging price while accounting for charging efficiency loss. Therefore, the reported cost provides a fair comparison even when the short- or medium-horizon case ends with a terminal SoC deficit. For the 24 h horizon, the original terminal SoC constraint is retained.

As shown in Table 6, the prediction horizon affects both the operating cost and computational burden. Shorter horizons reduce the optimization size and therefore require less decision time, but they provide less look-ahead information for future PV variation, load demand, and electricity-price changes. Meanwhile, shorter horizons also have difficulty in properly scheduling the batteries to meet the end-of-day terminal SoC target. Therefore, the terminal SoC deficits of the 6 h and 12 h cases are converted into equivalent charging compensation costs and included in the reported costs to avoid underestimating the cost of short-horizon operation. The 24 h horizon achieves the lowest cost because it can better coordinate battery charging and discharging over the whole daily cycle while satisfying the terminal SoC constraint. Across all tested horizons, the proposed MPC maintains a lower cost than the baseline SOCP-MPC, indicating that the advantage brought by the MMDG-AKRR-based PV correction is not limited to a specific prediction-horizon setting. In addition, the power-flow constraints remain satisfied in all tested cases, showing that the proposed online EMS preserves grid-security feasibility under different prediction-horizon configurations.

5.4 Key Results for Other Microgrid Cases↩︎

To avoid over-specialization to a specific network configuration, the numerical tests described above are further conducted on 10-bus, 33-bus, and 69-bus distribution grids, as shown in Fig. 13, Fig. 14, and Fig. 15.

The corresponding results are summarized in Table 7. Consistent conclusions can be drawn from these additional cases: (1) MPC-based methods generally achieve better economic performance compared with day-ahead EMS approaches due to their enhanced adaptability to real-time uncertainties; (2) Power-flow-constrained methods improve grid-security performance by explicitly accounting for branch-current and voltage constraints; (3) Cone-relaxed power-flow formulations are essential for enabling the EMS optimization to efficiently obtain high-quality solutions, whereas the nonlinear power-flow-constrained formulation incurs a significantly higher computational burden; and (4) The proposed MMDG-AKRR-assisted MPC-EMS delivers favorable overall performance by reducing operating costs while preserving grid-security performance and maintaining computational efficiency across the tested benchmark systems.

4pt

The tightness analysis for the cone-relaxation model can be found in 10.

6 Conclusion↩︎

A multi-modal dictionary-guided anisotropic kernel ridge regression (MMDG-AKRR) forecasting model-based MPC-EMS is proposed to address solar-forecasting deviations caused by uncertain weather and associated power-flow-constraint violations. By combining multiple representative base models through an online similarity-driven weighting mechanism, the proposed forecasting framework improves PV trajectory prediction under the represented multi-modal operating range. To ensure grid-secure operation, cone-relaxed power-flow constraints are explicitly embedded into the MPC-EMS formulation, enabling reliable and computationally efficient scheduling decisions. Numerical results demonstrate that the proposed MMDG-AKRR model achieves lower forecasting errors than the considered neural-network-based and kernel-based baselines under the tested regular, irregular, and multiple irregular PV profiles, making it suitable for online MPC applications that rely on frequent forecast updates. Extensive case studies on 10-bus, 18-bus, 33-bus, and 69-bus benchmark microgrids verify that the proposed MMDG-AKRR-assisted MPC-EMS achieves a favorable balance among operating cost, power-flow security, and computational efficiency. These results highlight the forecasting adaptability, operational effectiveness, and scalability-related computational tractability of the proposed framework across the tested benchmark systems for microgrid energy management under high renewable uncertainty.

7 Modal Dictionary Construction for MMDG-AKRR↩︎

The MMDG-AKRR framework requires a set of representative PV modal profiles to construct the base-model dictionary. In general, the modal dictionary can be obtained by either a data-driven clustering strategy or a typical-day scaling strategy.

For the clustering-based strategy, normalized daily PV profiles are grouped into \(K\) clusters using profile-distance-based methods such as K-means clustering. The samples assigned to the \(k\)th cluster are used to train the corresponding base AKRR model, and the cluster-average profile is used as the modal center curve for online similarity evaluation. The number of modals \(K\) determines the number of base models and is selected as a tradeoff between representation accuracy and computational complexity, which can be guided by validation performance.

For the typical-day scaling strategy, a representative clear-day profile \(\mathbf{p}_{\mathrm{typ}}\) is first selected from historical PV profiles, for example by identifying a group of highly similar clear-day profiles and taking their average. Then, \(K\) modal center profiles are generated by amplitude scaling: \[\begin{align} \bar{\mathbf{P}}_k = a_k \mathbf{P}_{\mathrm{typ}}, \quad k=1,\dots,K, \end{align}\] where \(a_k\) denotes the scaling factor used to cover different irradiance attenuation levels. Training samples around each scaled modal profile are then used to train the corresponding base AKRR model.

In the numerical study of this paper, the typical-day scaling strategy is adopted. This choice is suitable for the considered MPC-EMS application because PV generation has strong daily time-window regularity, and the equivalent PV energy over the prediction horizon is more influential for battery scheduling and grid-power exchange than local high-frequency profile fluctuations. Therefore, amplitude-scaled typical-day profiles provide a compact and practical modal dictionary for representing different PV generation levels in the tested EMS scenarios.

8 Multi-Start Multi-Step Validation Procedure↩︎

The objective function of the developed multi-start multi-step (MSMS) cross-validation is given in 35 . The multi-step property accounts for accumulated forecasting errors, while the multi-start property reduces sensitivity to initial points, which is essential for the MPC algorithm. \[\begin{align} \label{eq95MSMS} \min_{\sigma,\sigma_t,\lambda}\;\! \frac{1}{N_{v}}\!\sum_{v=1}^{N_v}\Big(\!\frac{1}{N_{s}}\!\sum_{s\in\Omega_s}\big(\!\frac{1}{N_{k}}\!\sum_{k=1}^{N_k} (\frac{x_k-(x_0+\sum_{i=1}^{k}f({\mathbf{x}}_n))}{x_k})^2\big)\!\Big) \end{align}\tag{35}\] where \(\mathbf{x}_n\) denotes the input sequence, \(x_0\) is the initial step, \(x_k\) represents the \(k\)th step sample. \(N_k\) is the number of forecasting steps. \(N_s\) is the number of initial points, belonging to the set \(\Omega_s\) and evenly distributed over a day. \(N_v\) denotes the number of the validation fold. For each fold, the validation set is divided into a training fold and a test fold. The hyperparameters of \(f\) are tuned on the training fold, while \(x_k\) and \(x_0\) are drawn from the test fold. The genetic algorithm (GA) is used as the solver for this validation process.

9 Quantitative Grid-Safety Statistics↩︎

To further quantify the grid-security performance shown in Fig. 10, the overlimit duration and maximum violation magnitude of the displayed bus voltage and branch current are reported in Table 8. The voltage violation is computed with respect to the lower and upper voltage limits, while the current violation is computed with respect to the allowable current magnitude of the displayed branch.

As shown in Table 8, the LP-based EMS cases lead to relatively large voltage and current violations, especially in terms of the maximum violation magnitudes. The SOCP-constrained day-ahead EMS reduces the maximum voltage and current violation magnitudes compared with the LP-based day-ahead EMS, although it still exhibits non-negligible overlimit durations under the tested dynamic simulation. In comparison, the online SOCP-MPC cases further suppress both voltage and current violations. In particular, the proposed MPC achieves the smallest voltage-violation duration and magnitude, with only \(0.88\%\) voltage-violation time and \(0.0035\) p.u. maximum voltage violation. It also eliminates the displayed branch-current overlimit, with \(0.00\%\) current-violation time and \(0.00\%\) maximum relative current violation. These quantitative results confirm that the proposed MPC effectively suppresses sustained and severe grid-security violations under the tested scenario.

10 Tightness Analysis of Cone-Relaxation↩︎

The tightness of the SOCP relaxation is further examined under the rolling MPC framework. The obtained gap record forms a two-dimensional matrix, where each row corresponds to one rolling MPC decision and each column corresponds to one look-ahead step. This allows both the worst-case relaxation gap and the rolling stability of the relaxation tightness to be evaluated.

To avoid inflation of the relative gap caused by numerically inactive branches, branch records with near-zero relaxation terms are screened out. Specifically, a branch record is regarded as active only when its denominator is no smaller than \(0.1\%\) of the maximum denominator among all branches at the same rolling decision and look-ahead step.

The extended tightness statistics are summarized in Table 9. In addition to the mean and standard deviation, the table reports the 95th percentile, the maximum gap, and the time, look-ahead step, and branch location of the worst-case gap. The maximum first-step gap is also reported because only the first step of each MPC optimization is physically applied to the microgrid. To examine whether relaxation-tightening oscillations appear during rolling optimization, a first-step oscillation index is computed as the maximum adjacent variation of the first-step relaxation gap over consecutive MPC decisions.

4pt

As shown in Table 9, the relaxation gaps remain small for most operating points. The 95th-percentile gap is \(3.62\%\), \(0.020\%\), \(0.025\%\), and \(0.026\%\) for the 10-, 18-, 33-, and 69-bus cases, respectively. The maximum full-horizon gap is \(3.98\%\), \(4.07\%\), \(3.49\%\), and \(2.46\%\), respectively. For the 18-, 33-, and 69-bus cases, the worst full-horizon gaps occur at relatively far look-ahead steps, while the actually implemented first-step gaps remain much smaller. For the 10-bus case, the worst gap occurs at the second look-ahead step and is close to the maximum first-step gap, but its magnitude remains below \(4\%\).

More importantly, the maximum first-step gap remains \(3.96\%\), \(0.236\%\), \(0.222\%\), and \(0.122\%\) for the four cases, respectively. Since only the first step of each rolling MPC solution is physically applied, these results indicate that the implemented dispatch decisions are only weakly affected by the cone relaxation. The first-step oscillation index is also limited to \(1.05\%\), \(0.225\%\), \(0.166\%\), and \(0.072\%\), respectively, showing that no evident relaxation-tightening oscillation appears during the rolling MPC process. These results support the satisfactory tightness of the SOCP relaxation for the proposed online EMS under the tested grid-connected scenarios, including the larger 69-bus system.

References↩︎