1 Introduction

With the rapid development of modern industry, the energy demands are increasing greatly. Compared with the traditional fossil fuels, the renewable energies are assumed to be clean, renewable and easy to be exploited. Considering the environment issues and resource crisis, the development of renewable energies has become an inevitable choice for countries all over the world [1]. Nowadays, the renewable energy generation industries, represented by wind energy and solar energy, have attracted considerable attention [2].

The active power generation may cause the frequency disturbances and changes of the grid due to the random and intermittent characteristics of renewable energies, and there are still great challenges in large-scale utilization of renewable energy. Therefore, smoothing the output power of renewable energies generation with uncertain characteristics is the basis to make efficient use of these energies.

For solving the output power fluctuation issue, energy storage equipment has been adopted to suppress the power fluctuation. In [3], a flywheel energy storage model based on permanent magnetic machines was proposed, and the positive impact on stability of power system was discussed. In [4], the energy storage system technologies were analyzed in the grid application. Restricted by the technology of manufacturing, however, the capacity and the lifetime of storage equipment are limited. When the fluctuations occur frequently or the output of renewable energies increases rapidly, the capacity of storage equipment will be insufficient, and meanwhile, the maintenance and replacement of energy storage equipment will become a great burden. Therefore, the adjustable loads are considered to solve the fluctuation problem.

In the past few years, there has been a growing interest in stochastic stability analysis. For instance, in [5], a stochastic small-signal stability analysis method is proposed for power system with wind power, and the impact of stochastic excitation is considered. In [6], the stochastic small-signal stability problem is discussed for the grid-connected photovoltaic system, and the influence of solar irradiation on system stability is also analyzed. In [7], a new concept is defined to study the small-signal stability of bulk power systems with wind generation. In [8], a probabilistic method is proposed to evaluate the possible effect of random wind power on the small-signal stability of power system. However, the stochastic characteristics of wind and solar irradiations are studied separately. Therefore, how to model the power system containing a variety of renewable energies which are coupled in the framework of stochastic system is an interesting and important research problem that has not been solved, which is one motivation of our research.

On the other hand, Markov model is widely utilized for describing the stochastic variation between finite states, and these variations can be input variables, state variables, or system parameters [9]. Therefore, based on Markov model, the analysis of the grid with stochastic renewable energy sources can no longer be limited to a single type of variable. In [10], a Markov process is used to describe the fluctuations caused by the local load, and then, the stochastic stability is analyzed, while the disturbances caused by power generation were not considered. In [11], a Markov process is adopted to describe the stochastically changing operating points. However, the above-mentioned results mainly concentrate on a single type of power generation. To realize the comprehensive utilization of the renewable energies, the multi-energy generation system has been studied [12], which contains a variety of renewable energies, and the stochastic characteristics of various renewable energies increase the complexity of the system. As mentioned before, Markov model is assumed to be a suitable tool for describing the stochastic fluctuations of renewable energies, especially for describing the system with a variety of renewable energies. To the best of the authors’ knowledge, there are few results considering this related topic; therefore, it is necessary and valuable to model and analyze multi-energy generation system based on Markov model.

Based on the above discussions, the main contributions of this paper are given as follows: To solve the renewable energies generation fluctuation issue, an adjustable loads control method is proposed by adjusting the power of adjustable loads reasonably. The multi-energy generation system with the integration of wind power and photovoltaic (PV) power is modeled in the framework of Markov theory. Furthermore, a new stochastic analysis method based on Markov model is given to illustrate the effectiveness of the proposed adjustable loads control strategy. To be noted that, this strategy not only can improve the utilization rate of renewable energies, but also can reduce the phenomenon of abandoning wind and solar energy.

The remaining of this paper is organized as follows. In Sect. 2, the multi-energy generation system with wind and PV is analyzed and Markov modeled. In Sect. 3, the numerical characteristics of wind speed and solar irradiation are analyzed, and the adjustable loads control strategy is designed according to the analysis results, then the stochastic stability of the multi-energy system is discussed based on stochastic Lyapunov theory. In Sect. 4, the effectiveness of the control strategy is verified by simulation results. Finally, conclusions are drawn.

2 Markov model of multi-energy generation system

To achieve the comprehensive utilization of renewable energies, the multi-energy generation system is proposed [13, 14], and the structure is flexible and unfixed, and it can be determined according to the local load demand and resource situation. In this paper, the typical photovoltaic and wind power are taken into account, and the adjustable load is also included [15,16,17]. In the following, the photovoltaic and wind power systems will be described in detail, respectively.

2.1 Photovoltaic system

A two-stage photovoltaic (PV) system mainly composes of four parts, a PV array as the dc source, a dc–dc converter with maximum power point tracking (MPPT) controller [18,19,20], a dc–ac inverter with a controller and the grid. The basic device of PV system is the PV cell that can directly convert solar energy into electricity. It follows from [21, 22] that the relationship between the current, the voltage, temperature and solar irradiance is described by

$$I = I_{pv} - I_{0} \left( {{\text{e}}^{{\frac{qV}{akT}}} - 1} \right) - \frac{{V + R_{\text{s}} I}}{{R_{\text{p}} }}$$
(1)

where Ipv and I0 are the solar irradiation current and leakage current, respectively. q is the electron charge, k is the Boltzmann constant, T is the temperature, \(\alpha\) is the ideality factor, and Rs and Rp are the parallel and series resistances.

In fact, many PV cells are usually combined into PV array when they are connected to the grid or in the practical application. The conversion efficiency of the PV cell is only around 21% [23]; therefore, some useful MPPT algorithms have been proposed [19, 24]. Usually, most of MPPT algorithms rely on dc–dc circuit. As the solar changes, the maximum power point moves, the corresponding duty ratio changes, and the PV output power, voltage and current also change. The relationship between input variables and output variables can be expressed as

$$\left\{ \begin{aligned} L\dot{i}_{\text{mp}} = v_{\text{mp}} - (1 - d(t))v_{\text{DC}} \hfill \\ C\dot{v}_{\text{DC}} = (1 - d(t))i_{\text{mp}} - i_{\text{DC}} \hfill \\ \end{aligned} \right.$$
(2)

where C and L are the capacitance and the inductance of the dc–dc circuit, respectively. d(t) is the duty ratio and \(v_{\text{mp}}\) and \(i_{\text{mp}}\) are the output voltage and current of the PV array, respectively. \(v_{\text{DC}}\) and \(i_{\text{DC}}\) are the voltage and current after dc–dc circuit, respectively.

PV array is dc power, to achieve grid connected, and it is necessary to convert dc into ac with the same voltage level as the grid-connected point. The three-phase voltage and current are transformed into the dq coordinate system. After the filter inductance \(L_{s}\), the voltages and currents are shown by following equations

$$\left\{ {\begin{array}{*{20}l} {L_{s} \dot{i}_{d} = e_{d} - Ri_{d} + \omega L_{s} i_{q} - v_{d} } \hfill \\ {L_{s} \dot{i}_{q} = e_{q} - Ri_{q} - \omega L_{s} i_{d} - v_{q} } \hfill \\ \end{array} } \right.$$
(3)

where id, iq, \(v_{d}\) and \(v_{q}\) are the currents and voltages in dq axis, respectively. ω is the voltage angular speed tracked by phase locked loop. ed and eq are the grid voltage in dq axis.

Without loss of generality, the PI controllers are applied to guarantee the photovoltaic system output voltage stability, and the current and voltage controllers are expressed by

$$\left\{ {\begin{array}{*{20}l} {\dot{x}_{1} = v_{{{\text{DC\_ref}}}} - v_{\text{DC}} } \hfill \\ {i_{{d{\text{\_ref}}}} = k_{1} (v_{{{\text{DC\_ref}}}} - v_{\text{DC}} ) + k_{2} x_{1} } \hfill \\ {\dot{x}_{2} = i_{{d{\text{\_ref}}}}^{s} - i_{d} } \hfill \\ {u_{d} = k_{3} (i_{{d{\text{\_ref}}}}^{s} - i_{d} ) + k_{4} x_{2} - \omega L_{s} i_{q} + e_{d} } \hfill \\ {\dot{x}_{3} = i_{{q{\text{\_ref}}}}^{s} - i_{q} } \hfill \\ {u_{q} = k_{5} (i_{{q{\text{\_ref}}}}^{s} - i_{q} ) + k_{6} x_{3} + \omega L_{s} i_{d} + e_{q} } \hfill \\ \end{array} } \right.$$
(4)

where \(v_{{{\text{DC\_ref}}}}\) is the reference voltage of the capacitor and \(i_{{d{\text{\_ref}}}}^{s}\) and \(i_{{q{\text{\_ref}}}}^{s}\) are the reference current in d axis and q axis, respectively. k1, k3 and k5 are the proportionality coefficients of the controller. k2, k4 and k6 are the integral coefficients of the controller. x1, x2 and x3 are intermediate variables.

According to formula (2)–(4), the closed-loop system in state space could be expressed by

$$\left[ {\begin{array}{*{20}l} {\dot{i}_{\text{mp}} } \hfill \\ {\dot{v}_{\text{DC}} } \hfill \\ {i_{d} } \hfill \\ {\dot{i}_{q} } \hfill \\ {\dot{x}_{1} } \hfill \\ {\dot{x}_{2} } \hfill \\ {\dot{x}_{3} } \hfill \\ \end{array} } \right] = \left[ {\begin{array}{*{20}l} 0 \hfill & { - \frac{d(t)}{L}} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ {\frac{1 - d(t)}{C}} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & { - \frac{{k_{1} k_{3} }}{{L_{s} }}} \hfill & { - \frac{{k_{3} }}{{L_{s} }}} \hfill & 0 \hfill & {\frac{{k_{2} k_{3} }}{{L_{s} }}} \hfill & {\frac{{k_{4} }}{{L_{s} }}} \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \frac{{k_{5} }}{{L_{s} }}} \hfill & 0 \hfill & 0 \hfill & {\frac{{k_{6} }}{{L_{s} }}} \hfill \\ 0 \hfill & { - 1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & { - k_{1} } \hfill & { - 1} \hfill & 0 \hfill & {k_{2} } \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - 1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ \end{array} } \right] \times \left[ {\begin{array}{*{20}l} {i_{\text{mp}} } \hfill \\ {v_{\text{DC}} } \hfill \\ {i_{d} } \hfill \\ {i_{q} } \hfill \\ {x_{1} } \hfill \\ {x_{2} } \hfill \\ {x_{3} } \hfill \\ \end{array} } \right] + \left[ {\begin{array}{*{20}l} 0 \hfill & 0 \hfill & {\frac{1}{L}} \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \frac{1}{C}} \hfill & 0 \hfill & 0 \hfill \\ {\frac{1}{{L_{s} }}} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & {\frac{{k_{1} k_{3} }}{{L_{s} }}} \hfill & 0 \hfill \\ 0 \hfill & {\frac{1}{{L_{s} }}} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & {\frac{{k_{5} }}{{L_{s} }}} \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 1 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & {k_{1} } \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 1 \hfill \\ \end{array} } \right] \times \left[ {\begin{array}{*{20}l} {e_{d} } \hfill \\ {e_{q} } \hfill \\ {v_{\text{mp}} } \hfill \\ {i_{\text{DC}} } \hfill \\ {v_{{{\text{DC\_ref}}}} } \hfill \\ {i_{{{\text{q\_ref}}}} } \hfill \\ \end{array} } \right]$$
(5)

2.2 Wind power system

A grid-connected wind power generation system mainly consists of wind wheel, gear box, generator, converter and the grid [25,26,27]. As the wind speed changes, many variables also change, such as the torque and the output power. The torque obtained from the wind turbine is given by [28]

$$T_{e} = \frac{{0.5\pi \rho R^{2} V^{3} C_{p} }}{{\omega_{t} }}$$
(6)

where Te represents the torque, \(\rho\) represents the air density, R represents the radius of the wind wheel, V represents the wind speed, and \(\omega_{t}\) represents the angular speed of the wind turbine. \(C_{p}\) is expressed as

$$\left\{ {\begin{array}{*{20}l} {C_{p} = 0.22\left( {\frac{116}{\lambda } - 0.4\beta - 5} \right){\text{e}}^{{{\raise0.7ex\hbox{${ - 12.5}$} \!\mathord{\left/ {\vphantom {{ - 12.5} {\lambda_{i} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${\lambda_{i} }$}}}} } \hfill \\ {\frac{1}{{\lambda_{i} }} = \frac{1}{\lambda + 0.08\beta } + \frac{0.035}{{\beta^{3} + 1}}} \hfill \\ \end{array} } \right.$$
(7)

where \(\beta\) represents the pitch angle, and \(\lambda\) represents the tip speed ratio.

In [29,30,31], the stator flux linkage model of induction generator is presented as

$$\left\{ {\begin{array}{*{20}l} {v_{{d{\text{s}}}} = R_{\text{s}} i_{{d{\text{s}}}} - \omega_{\text{e}} \psi_{{q{\text{s}}}} + \dot{\psi }_{{d{\text{s}}}} } \hfill \\ {v_{{q{\text{s}}}} = R_{\text{s}} i_{{q{\text{s}}}} + \omega_{\text{e}} \psi_{{d{\text{s}}}} + \dot{\psi }_{{q{\text{s}}}} } \hfill \\ {\psi_{{d{\text{s}}}} = L_{\text{s}} i_{{d{\text{s}}}} + L_{\text{m}} i_{{d{\text{r}}}} } \hfill \\ {\psi_{{q{\text{s}}}} = L_{\text{s}} i_{{q{\text{s}}}} + L_{\text{m}} i_{{q{\text{r}}}} } \hfill \\ \end{array} } \right.$$
(8)

The rotor flux linkage model of induction generator is given by

$$\left\{ {\begin{array}{*{20}l} {v_{{d{\text{r}}}} = R_{\text{r}} i_{{d{\text{r}}}} - (\omega_{\text{e}} - \omega_{\text{r}} )\psi_{{q{\text{r}}}} + \dot{\psi }_{{d{\text{r}}}} } \hfill \\ {v_{{q{\text{r}}}} = R_{\text{r}} i_{{q{\text{r}}}} + (\omega_{\text{e}} - \omega_{\text{r}} )\psi_{{d{\text{r}}}} + \dot{\psi }_{{q{\text{r}}}} } \hfill \\ {\psi_{{d{\text{r}}}} = L_{\text{r}} i_{{d{\text{r}}}} + L_{\text{m}} i_{{d{\text{s}}}} } \hfill \\ {\psi_{{q{\text{r}}}} = L_{\text{r}} i_{{q{\text{r}}}} + L_{\text{m}} i_{{q{\text{s}}}} } \hfill \\ \end{array} } \right.$$
(9)

where \(v_{{d{\text{s}}}}\), \(v_{{q{\text{s}}}}\), \(v_{{d{\text{r}}}}\) and \(v_{{q{\text{r}}}}\) are the components of the stator and rotor terminal voltage on the d axis and the q axis, respectively. \(\psi_{{d{\text{s}}}}\), \(\psi_{{q{\text{s}}}}\), \(\psi_{{d{\text{r}}}}\) and \(\psi_{{q{\text{r}}}}\) are the components of stator and rotor flux linkage on the d axis and the q axis, respectively. \(i_{{d{\text{r}}}}\), \(i_{{q{\text{r}}}}\), \(i_{{d{\text{s}}}}\) and \(i_{{q{\text{s}}}}\) are the components of rotor and stator current on the d axis and the q axis, respectively. Rs and Rr are the stator and rotor resistances, respectively. \(\omega_{\text{e}}\) and \(\omega_{\text{r}}\) represent the synchronous speed and rotor speed, respectively.

Similarly, the PI controllers are applied in the wind power system, and the current controllers are expressed by

$$\left\{ {\begin{array}{*{20}l} {\dot{x}_{7} = i_{{d{\text{r\_ref}}}}^{w} - i_{{d{\text{r}}}} } \hfill \\ {\dot{x}_{8} = i_{{q{\text{r\_ref}}}}^{w} - i_{{q{\text{r}}}} } \hfill \\ {v_{{d{\text{r}}}} = k_{7} (i_{{d{\text{r\_ref}}}}^{w} - i_{{d{\text{r}}}} ) + k_{8} x_{7} + R_{\text{r}} i_{{d{\text{r}}}} - \omega_{\text{e}} \psi_{{q{\text{r}}}} } \hfill \\ {v_{{q{\text{r}}}} = k_{9} (i_{{q{\text{r\_ref}}}}^{w} - i_{{q{\text{r}}}} ) + k_{10} x_{8} + R_{\text{r}} i_{{q{\text{r}}}} + \omega_{\text{e}} \psi_{{d{\text{r}}}} } \hfill \\ \end{array} } \right.$$
(10)

where \(i_{{d{\text{r\_ref}}}}^{w}\) is the reference current in d axis and \(i_{{q{\text{r\_ref}}}}^{w}\) is the reference current in q axis. k7 and k9 are the proportionality coefficients, k8 and k10 are the integral coefficients of the controller, and x7 and x8 are intermediate variables.

According to formula (8)–(10), the closed-loop wind power system could be expressed by

$$\begin{aligned} \left[ {\begin{array}{*{20}l} {\dot{\psi }_{{d{\text{s}}}} } \hfill \\ {\dot{\psi }_{{q{\text{s}}}} } \hfill \\ {\dot{i}_{{d{\text{r}}}} } \hfill \\ {\dot{i}_{{q{\text{r}}}} } \hfill \\ {\dot{x}_{7} } \hfill \\ {\dot{x}_{8} } \hfill \\ \end{array} } \right] & = \left[ {\begin{array}{*{20}l} { - \,\frac{{R_{\text{s}} }}{{L_{\text{s}} }}} \hfill & {\omega_{\text{e}} } \hfill & { - \,R_{\text{s}} } \hfill & {\frac{{L_{\text{m}} R_{\text{s}} }}{{L_{\text{s}} }}} \hfill & 0 \hfill & 0 \hfill \\ { - \,\omega_{\text{e}} } \hfill & { - \frac{{R_{\text{s}} }}{{L_{\text{s}} }}} \hfill & 0 \hfill & { - \,R_{\text{s}} } \hfill & {\frac{{L_{\text{m}} R_{\text{s}} }}{{L_{\text{s}} }}} \hfill & 0 \hfill \\ {\frac{{R_{\text{s}} L_{\text{m}} }}{{L_{\text{s}} (L_{\text{s}} L_{\text{r}} - L_{\text{m}}^{2} )}}} \hfill & { - \,\frac{{(\omega_{\text{e}} + \omega_{\text{r}} )L_{\text{m}} }}{{L_{\text{s}} L_{\text{r}} - L_{\text{m}}^{2} }}} \hfill & { - \frac{{R_{\text{s}} L_{\text{m}}^{2} + k_{7} L_{\text{s}}^{2} }}{{L_{\text{s}} (L_{\text{s}} L_{\text{r}} - L_{\text{m}}^{2} )}}} \hfill & 0 \hfill & {\frac{{L_{\text{s}} k_{8} }}{{L_{\text{s}} L_{\text{r}} - L_{\text{m}}^{2} }}} \hfill & 0 \hfill \\ {\frac{{(\omega_{\text{e}} + \omega_{\text{r}} )L_{\text{m}} }}{{L_{\text{s}} L_{\text{r}} - L_{\text{m}}^{2} }}} \hfill & {\frac{{R_{\text{s}} L_{\text{m}} }}{{L_{\text{s}} (L_{\text{s}} L_{\text{r}} - L_{\text{m}}^{2} )}}} \hfill & 0 \hfill & { - \,\frac{{R_{\text{s}} L_{\text{m}}^{2} + k_{9} L_{\text{s}}^{2} }}{{L_{\text{s}} (L_{\text{s}} L_{\text{r}} - L_{\text{m}}^{2} )}}} \hfill & 0 \hfill & {\frac{{L_{\text{s}} k_{10} }}{{L_{\text{s}} L_{\text{r}} - L_{\text{m}}^{2} }}} \hfill \\ 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill \\ \end{array} } \right] \times \left[ {\begin{array}{*{20}l} {\psi_{{d{\text{s}}}} } \hfill \\ {\psi_{{q{\text{s}}}} } \hfill \\ {i_{{d{\text{r}}}} } \hfill \\ {i_{{q{\text{r}}}} } \hfill \\ {x_{7} } \hfill \\ {x_{8} } \hfill \\ \end{array} } \right] \\ & \quad + \left[ {\begin{array}{*{20}l} 1 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & 1 \hfill & 0 \hfill & 0 \hfill \\ { - \,\frac{{L_{\text{m}} }}{{L_{\text{s}} L_{\text{r}} - L_{\text{m}}^{2} }}} \hfill & 0 \hfill & {\frac{{L_{\text{s}} k_{7} }}{{L_{\text{s}} L_{\text{r}} - L_{\text{m}}^{2} }}} \hfill & 0 \hfill \\ 0 \hfill & { - \,\frac{{L_{\text{m}} }}{{L_{\text{s}} L_{\text{r}} - L_{\text{m}}^{2} }}} \hfill & 0 \hfill & {\frac{{L_{\text{s}} k_{9} }}{{L_{\text{s}} L_{\text{r}} - L_{\text{m}}^{2} }}} \hfill \\ 0 \hfill & 0 \hfill & 1 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & 1 \hfill \\ \end{array} } \right] \times \left[ {\begin{array}{*{20}l} {e_{d} } \hfill \\ {e_{q} } \hfill \\ {i_{{d{\text{r\_ref}}}} } \hfill \\ {i_{{q{\text{r\_ref}}}} } \hfill \\ \end{array} } \right] \\ \end{aligned}$$
(11)

2.3 Markov model

In general, a Markov model consists of a Markov chain, a probability transition matrix and a state space equation [32, 33]. A Markov chain with finite number of modes is defined to describe the variation of solar irradiation and wind speed, and the solar data and wind data are clustered, then modes of solar irradiation and modes of wind speed are obtained, and every mode includes a cluster center [34].

To analyze the multi-energy generation system, the coupled Markov model of solar and wind are considered. Let there be m solar irradiation modes and n wind speed modes, the multi-energy generation system would have \(m \times n\) modes. The relationship between multi-energy generation system modes, wind modes and solar modes is shown in Table 1.

Table 1 Relationship between the multi-energy generation system modes, wind modes and solar modes
Full size table

Based on the clustering results, the wind speed, the solar irradiation, the related duty ratio d(r(t)) of PV and the rotor speed of ωr wind power in each mode are calculated. Substituting the duty ratio and rotor speed values of each mode into the formula (5) and (11), the continuous linear Markov model of multi-energy generation system is obtained which is as follows:

$$\dot{x}(t) = A(r(t))x(t) + B(r(t))u(t)$$
(12)

where \(x = [\begin{array}{*{20}c} {i_{\text{mp}} } & {v_{\text{DC}} } & {i_{d} } & {i_{q} } & {x_{1} } & {x_{2} } & {x_{3} } & {\psi_{{d{\text{s}}}} } & {\psi_{dq} } & {i_{{d{\text{r}}}} } & {i_{{q{\text{r}}}} } & {x_{7} } & {x_{8} } \\ \end{array} ]\) and \(u = [\begin{array}{*{20}l} {v_{\text{mp}} } \hfill & {i_{\text{DC}} } \hfill & {v_{{{\text{DC\_ref}}}} } \hfill & {i_{{q{\text{\_ref}}}}^{s} } \hfill & {e_{d} } \hfill & {e_{q} } \hfill & {i_{{d{\text{r\_ref}}}} } \hfill & {i_{{q{\text{r\_ref}}}} } \hfill \\ \end{array} ]\) are the state vector and the input vector, respectively. \(r(t) = 1,2, \ldots ,m \times n\) represents different modes, and A(r(t)) and B(r(t)) are the matrices corresponding to the different modes.

By quantifying and clustering the data and using Bayesian estimation and inference [35, 36], the transition probability of solar and wind is calculated. Since the wind speed and solar irradiation are independent, the probability that they both occur at the same time are the product of their respective probability

$$p_{mn} = p_{{_{ij} }}^{s} \times p_{{_{kh} }}^{w}$$
(13)

where pmn is the transition probability of the multi-energy generation system and \(p_{{_{ij} }}^{s}\) and \(p_{{_{kh} }}^{w}\) are the transition probability of solar irradiation and wind speed, respectively.

According to the transition probability [37], the transition rate matrix \(\varPi\) of the multi-energy generation system can be obtained. To be noted that the diagonal elements should be non-positive to satisfy the basic Markov theory rule that the sum of the elements of a row is zero.

$$\varPi {\text{ = \{ }}\pi_{mn} {\text{\} ,}}\pi_{mn} = \left\{ {\begin{array}{*{20}l} {P_{mn} ,} \hfill & {m \ne n} \hfill \\ {\sum\limits_{m \ne n} {P_{mn} - 1} ,} \hfill & {m = n} \hfill \\ \end{array} } \right.$$
(14)

Remark 1

In this paper, the coupled Markov model of the multi-energy generation system is proposed, and the developed coupled Markov model not only reflects the stochastic variables, such as wind speed and solar irradiation by finite modes and the probability of transition between each mode, but also can make more detailed and accurate description of the current and power change within the multi-energy generation system by introducing stochastic variables into the state space equation [38]. Based on the developed Markov model, the multi-energy generation system can be further analyzed in detail.

3 Adjustable loads control strategy and stochastic stability analysis

In this section, to weaken the power generation fluctuation and improve the utilization of the renewable energies, the numerical characteristics of the multi-energy generation system are analyzed and an adjustable control strategy is proposed to guarantee the continuous stability of the output power. Furthermore, the stochastic stability of multi-energy generation system is discussed based on Markov model.

3.1 The numerical characteristics of renewable energies

The power generation of the multi-energy generation system considered in this paper entirely depends on renewable energies. Commonly, energy storage device plays a smoothing role in power fluctuations, and the storage capacity is a significant parameter, and if the capacity setting is too small, the smoothing function is weak, and if the capacity setting is too large, the device cannot be fully utilized.

To obtain a reasonable capacity of energy storage device, the numerical characteristics of renewable energies should be analyzed. According to Markov theory, an irreducible aperiodic Markov chain exists at least one stationary distribution \(\pi = (\pi_{1} ,\pi_{2} , \ldots ,\pi_{s} )\) satisfying the following relations [39, 40]

$$\sum\limits_{i = 1}^{s} {\pi {}_{i}} = 1,\quad \pi_{i} \ge 0$$
(15)
$$\mathop {\text{Lim}}\limits_{n \to \infty } \left\| {P^{n} X^{n} - \pi } \right\| = 0$$
(16)

Based on the Markov model, the stationary distribution can be obtained, and the output power expectation of the multi-energy generation system can be calculated

$$\hat{P} = \sum\limits_{j = 1}^{j = s} {\pi_{j} y_{j} } + \sum\limits_{i = 1}^{i = s} {\pi_{i} x_{i} }$$
(17)

where \(\pi_{j}\) and \(\pi_{i}\) represent the probability of solar and wind, respectively. yi and xi represent the output power of PV and wind power in the corresponding mode, respectively.

The expectation value is considered as the system stability output, and the standard deviation reflects the extent of fluctuations, and it can be calculated by

$$d = \sqrt {\sum\limits_{i = 1,j = 1}^{i = s,j = s} {\pi_{i} \pi_{j} (x_{i} + y_{i} - \hat{P})^{2} } }$$
(18)

Remark 2

The capacity of energy storage equipment is mainly dependent on the power fluctuation. Based on the developed Markov model, the numerical characteristics are obtained. According to the numerical characteristics of the renewable energies, the fluctuation of renewable energies power generation can be described. Standard deviation is a suitable characterization index. Therefore, the capacity value of the energy storage device is set more than the standard deviation value.

3.2 The adjustable loads control strategy

Constrained by the capacity and storage conditions, the regulation ability of storage device is low, especially when the output power remains high value or the output power fluctuates continuously. A better scheme is to take advantage of the adjustable loads to consume the excess energy. Hydrogen and hot water are easy to store and transport [41]; hence, the extra electric energy can be used to produce hydrogen and heating, according to the demands. Nowadays, the hydrogen production is mainly manufactured by electrolysis of sodium hydroxide solutions [42]. The electrochemical reactions at the cathode and anode are given as follows

$${\text{Cathode:}}\quad 2{\text{H}}_{ 2} {\text{O}} + 4e_{ - } \to {\text{H}}_{2} \uparrow + 4{\text{OH}}_{ - }$$
(19)
$${\text{Anode}}:\quad 4{\text{OH}}_{ - } \to {\text{O}}_{2} \uparrow + 2{\text{H}}_{ 2} {\text{O}} + 4e_{ - }$$
(20)

The relationship between the hydrogen production amount and the electricity power consumption is described by

$$M_{{ - {\text{H}}_{ 2} }} { = }1.05 \times 10^{ - 8} \left( {\frac{{P_{{ - {\text{H}}_{ 2} }} }}{{v_{\text{DC}} }}} \right)$$
(21)

where \(M_{{ - {\text{H}}_{ 2} }}\) represents the quality of hydrogen and \(P_{{ - {\text{H}}_{ 2} }}\) is the electric power consumed by hydrogen production.

Ignoring the heat dissipation, the relationship between heating water and electricity power consumption can be formulated as

$$M_{\text{w}} C(T_{\text{tar}} - T_{\text{ref}} ) = \int_{{t_{1} }}^{{t_{2} }} {P(t)} {\text{d}}t$$
(22)

where Mw represents the quality of water, C is the specific heat of water, \(T_{\text{ref}}\) is the room temperature, \(T_{\text{tar}}\) is the target temperature, and P(t) is the consumed power.

Considering the stochastic fluctuation of output power Pgen caused by renewable energies, to guarantee the continuous and constant power transmission PTra to the grid, the following controller is designed for the flexible loads

$$P_{\text{con}} = (k_{p} + k_{i} /s)(i_{\text{tra}} - i_{\text{ref}} )$$
(23)

where Pcon is the amount of the flexible loads consumption, iref is the reference current transmitted to the grid, and itra is the actual current transmitted to the grid.

Considering the power regulation range and the ramp rate limit, the constraints of the flexible loads are formulated as

$$\left\{ {\begin{array}{*{20}l} {0 \le P_{\text{con}} \le P_{{{\text{con\_max}}}} } \hfill \\ {\lambda_{\text{lower}} \le \frac{{{\text{d}}P_{\text{con}} }}{{{\text{d}}t}} \le \lambda_{\text{upper}} } \hfill \\ \end{array} } \right.$$
(24)

where \(P_{{{\text{con\_max}}}}\) is the maximum power consumption and \(\lambda_{\text{up}}\) and \(- \,\lambda_{\text{down}}\) are the upper limit and lower limit of the ramp constraints, respectively.

When the constraints are satisfied, the consumption of adjustable loads depends on the power fluctuations of the renewable energies. According to (21)–(23), the relationship between the current and the production amounts of the hot water and the hydrogen is given by

$$\frac{{v_{\text{DC}} \times M_{{ - {\text{H}}_{2} }} }}{{1.05 \times 10^{ - 8} }} + M_{\text{w}} C\frac{{{\text{d}}T}}{{{\text{d}}t}} = \left( {k_{p} + \frac{{k_{i} }}{s}} \right)\left( {i_{\text{tra}} - i_{\text{ref}} } \right)$$
(25)

It is assumed that the reactive power is constant and the active power control is mainly considered in this paper, and the relationship between flexible loads current and the output current in d axis is as follows

$$i_{\text{dr}} = i_{\text{tra}} + i_{\text{loads}} = i_{\text{tra}} + \frac{{P_{\text{loads}} }}{{e_{d} }} = i_{\text{tra}} + \left( {k_{p} + \frac{{k_{i} }}{s}} \right)\frac{{\left( {i_{\text{tra}} - i_{\text{ref}} } \right)}}{{e_{d} }}$$
(26)

According to (12) and (25), the Markov model of the multi-energy system with adjustable loads control can be re-described as

$$\dot{x}(t) = A_{i} x(i) + B_{i} u(i)$$
(27)

where \(u(i) = [\begin{array}{*{20}c} {v_{\text{mp}} } & {i_{\text{DC}} } & {v_{{{\text{DC\_ref}}}} } & {i_{{q{\text{\_ref}}}}^{s} } & {e_{d} } & {e_{q} } & {i_{{d{\text{r\_ref}}}} } & {i_{{q{\text{r\_ref}}}} } \\ \end{array} ]\) and \(x(i) = [\begin{array}{*{20}c} {i_{\text{mp}} } & {v_{\text{DC}} } & {i_{d} } & {i_{q} } & {x_{1} } & {x_{2} } & {x_{3} } & {\psi_{{d{\text{s}}}} } & {\psi_{{d{\text{q}}}} } & {i_{{d{\text{r}}}} } & {i_{{q{\text{r}}}} } & {x_{7} } & {x_{8} } \\ \end{array} ]\) are the input vector and the state vector, respectively.

Remark 3

The adjustable loads play a significant role in restraining the fluctuation of renewable energy generation. The advantage is that it can fully use the wind and solar energy. To be pointed out that the manufacturing hydrogen and heating water are just examples in this paper, the actual configuration of adjustable loads should be based on local requirements, and the adjustable loads should be easy to adjust, and their products should be convenient for transportation and storage.

3.3 Stochastic stability analysis

Theorem 1

The Markov model of multi-energy generation system (27) with adjustable loads control is stochastically stable if and only if there exists a set of positive definite symmetric matrices\(P_{i} > 0,i,j \in S = \{ 1,2, \ldots m \times n\}\), satisfying

$$A_{i}^{T} P_{i} + P_{i} A_{i} + \sum\limits_{i = 1}^{s} {\pi_{ij} P_{j} } < 0$$
(28)

where Aiis the linearized state matrix at different operating modes,\(\pi_{ij}\)is the transition rate, and Piis a positive definite symmetric matrix.

Proof

Construct the following stochastic Lyapunov function

$$V(x(t),r(t) = i) = V(x(t),i) = x^{T} (t)P_{i} x(t)$$
(29)

For all \(x(t) \ne 0\), we have \(V(x(t),i) > 0\). The weak infinitesimal operator \(\ell V(x(t),i)\) of the Lyapunov function is obtained

$$\begin{aligned} \ell V\left( {x(t),i} \right) & = \mathop {\text{Lim}}\limits_{\delta \to 0} \frac{1}{\delta }\left\{ {E\left[ {V\left( {x\left( {\left( {t + \delta } \right),r(t + \delta )} \right)|x(t),i} \right)} \right] - V(x,i)} \right\} \\ & = \sum\limits_{j = 1}^{s} {\pi_{ij} V(x(r(t),r(t + \delta ))) + x^{T} (r(t))A_{i}^{T} \frac{\partial }{\partial x}} V(x(r(t),i)) \\ & = x^{T} \left[ {\sum\limits_{j = 1}^{s} {\pi_{ij} P_{j} + A_{i}^{T} P_{i} + P_{i} A_{i} } } \right]x \\ \end{aligned}$$
(30)

It follows from condition (27) that

$$\ell V(x(t),i) = x^{T} \left[ {\sum\limits_{j = 1}^{s} {\pi_{ij} P_{j} + A_{i}^{T} P_{i} + P_{i} A_{i} } } \right]x < 0$$
(31)

The multi-energy generation system (27) is stochastically stable if condition (28) is satisfied, and the detailed calculation is provided in the case study.

4 Case study

In this paper, the multi-energy generation system containing a wind power and a photovoltaic is considered. The type of photovoltaic cell is SPR-305-WHT, and a PV module consists of 96 photovoltaic cells. The 5 photovoltaic modules form a string, and the 66 strings form a photovoltaic array, and the specific parameters are shown in Table 2.

Table 2 Parameters of PV cell
Full size table

To make a reasonable coordination, the wind power and photovoltaic power generation are set as the same capacity level, and the specific parameters are provided in Table 3.

Table 3 Parameters of generator
Full size table

Four discrete states are defined in this paper by using a K-means clustering analysis of the collected data, the Markov modules for solar and wind speed are shown in Tables 4 and 5, respectively.

Table 4 Markov modules for solar irradiation
Full size table
Table 5 Markov modules for wind speed
Full size table

The coefficients of PI controllers and the positive definite symmetric matrix Pi are set as

$$\begin{aligned} k_{1} & = 0.75,k_{3} = k_{5} = 0.5,k_{2} = 10,k_{4} = k_{6} = 60 \\ k_{7} & = k_{9} = 0.155,k_{8} = k_{10} = 15,P_{i} = I_{13 \times 13} \\ \end{aligned}$$

The matrix Ai of the multi-energy system in sixteen modes is provided as follows

$$A_{i} = \left[ {\begin{array}{*{20}c} {\psi_{m} } & 0 \\ 0 & {\varphi_{n} } \\ \end{array} } \right],\quad \{ i = 1,2, \ldots ,16|m,n = 1,2, \ldots ,4\}$$
$$\psi_{1} = \left[ {\begin{array}{*{20}l} 0 \hfill & { - \,0.187} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ {0.072} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & { - \,1} \hfill & { - \,1.333} \hfill & 0 \hfill & {13.33} \hfill & {160} \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \,1.333} \hfill & 0 \hfill & 0 \hfill & {160} \hfill \\ 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & { - \,0.75} \hfill & { - \,1} \hfill & 0 \hfill & {10} \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ \end{array} } \right]$$
$$\psi_{2} = \left[ {\begin{array}{*{20}l} 0 \hfill & { - \,0.213} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ {0.067} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & { - \,1} \hfill & { - \,1.333} \hfill & 0 \hfill & {13.33} \hfill & {160} \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \,1.333} \hfill & 0 \hfill & 0 \hfill & {160} \hfill \\ 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & { - \,0.75} \hfill & { - \,1} \hfill & 0 \hfill & {10} \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ \end{array} } \right]$$
$$\psi_{3} = \left[ {\begin{array}{*{20}l} 0 \hfill & { - \,0.24} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ {0.061} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & { - \,1} \hfill & { - \,1.333} \hfill & 0 \hfill & {13.33} \hfill & {160} \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \,1.333} \hfill & 0 \hfill & 0 \hfill & {160} \hfill \\ 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & { - \,0.75} \hfill & { - \,1} \hfill & 0 \hfill & {10} \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ \end{array} } \right]$$
$$\psi_{4} = \left[ {\begin{array}{*{20}l} 0 \hfill & { - \,0.267} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ {0.056} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & { - \,1} \hfill & { - \,1.333} \hfill & 0 \hfill & {13.333} \hfill & {160} \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \,1.333} \hfill & 0 \hfill & 0 \hfill & {160} \hfill \\ 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & { - \,0.75} \hfill & { - \,1} \hfill & 0 \hfill & {10} \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ \end{array} } \right]$$
$$\varphi_{1} = \left[ {\begin{array}{*{20}l} { - \,0.00198} \hfill & 1 \hfill & { - \,0.00607} \hfill & {0.00573} \hfill & 0 \hfill & 0 \hfill \\ { - \,1} \hfill & { - \,0.00198} \hfill & 0 \hfill & { - \,0.00607} \hfill & {0.00573} \hfill & 0 \hfill \\ {0.00588} \hfill & { - \,5.042} \hfill & { - \,0.50523} \hfill & 0 \hfill & {47.2473} \hfill & 0 \hfill \\ {5.042} \hfill & {0.00588} \hfill & 0 \hfill & { - \,0.50523} \hfill & 0 \hfill & {47.2473} \hfill \\ 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill \\ \end{array} } \right]$$
$$\varphi_{2} = \left[ {\begin{array}{*{20}l} { - \,0.00198} \hfill & 1 \hfill & { - \,0.00607} \hfill & {0.00573} \hfill & 0 \hfill & 0 \hfill \\ { - \,1} \hfill & { - \,0.00198} \hfill & 0 \hfill & { - \,0.00607} \hfill & {0.00573} \hfill & 0 \hfill \\ {0.00588} \hfill & { - \,5.268} \hfill & { - \,0.50523} \hfill & 0 \hfill & {47.2473} \hfill & 0 \hfill \\ {5.268} \hfill & {0.00588} \hfill & 0 \hfill & { - \,0.50523} \hfill & 0 \hfill & {47.2473} \hfill \\ 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill \\ \end{array} } \right]$$
$$\varphi_{2} = \left[ {\begin{array}{*{20}l} { - \,0.00198} \hfill & 1 \hfill & { - \,0.00607} \hfill & {0.00573} \hfill & 0 \hfill & 0 \hfill \\ { - \,1} \hfill & { - \,0.00198} \hfill & 0 \hfill & { - \,0.00607} \hfill & {0.00573} \hfill & 0 \hfill \\ {0.00588} \hfill & { - \,5.268} \hfill & { - \,0.50523} \hfill & 0 \hfill & {47.2473} \hfill & 0 \hfill \\ {5.268} \hfill & {0.00588} \hfill & 0 \hfill & { - \,0.50523} \hfill & 0 \hfill & {47.2473} \hfill \\ 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill \\ \end{array} } \right]$$
$$\varphi_{3} = \left[ {\begin{array}{*{20}l} { - \,0.00198} \hfill & 1 \hfill & { - \,0.00607} \hfill & {0.00573} \hfill & 0 \hfill & 0 \hfill \\ { - \,1} \hfill & { - \,0.00198} \hfill & 0 \hfill & { - \,0.00607} \hfill & {0.00573} \hfill & 0 \hfill \\ {0.00588} \hfill & { - \,5.922} \hfill & { - \,0.50523} \hfill & 0 \hfill & {47.2473} \hfill & 0 \hfill \\ {5.922} \hfill & {0.00588} \hfill & 0 \hfill & { - \,0.50523} \hfill & 0 \hfill & {47.2473} \hfill \\ 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill \\ \end{array} } \right]$$
$$\varphi_{4} = \left[ {\begin{array}{*{20}l} { - \,0.00198} \hfill & 1 \hfill & { - \,0.00607} \hfill & {0.00573} \hfill & 0 \hfill & 0 \hfill \\ { - \,1} \hfill & { - \,0.00198} \hfill & 0 \hfill & { - \,0.00607} \hfill & {0.00573} \hfill & 0 \hfill \\ {0.00588} \hfill & { - \,6.506} \hfill & { - \,0.50523} \hfill & 0 \hfill & {47.2473} \hfill & 0 \hfill \\ {6.506} \hfill & {0.00588} \hfill & 0 \hfill & { - \,0.50523} \hfill & 0 \hfill & {47.2473} \hfill \\ 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & 0 \hfill & { - \,1} \hfill & 0 \hfill & 0 \hfill \\ \end{array} } \right]$$

As discussed before, the transition probabilities of solar irradiation and wind speed are obtained as follows:

$$\varPi^{s} = \left[ {\begin{array}{*{20}l} { - \,1} \hfill & {0.67} \hfill & {0.33} \hfill & 0 \hfill \\ {0.0294} \hfill & { - \,0.6471} \hfill & {0.3824} \hfill & {0.2353} \hfill \\ {0.0218} \hfill & {0.3478} \hfill & { - \,0.5435} \hfill & {0.1739} \hfill \\ {0.0588} \hfill & {0.2942} \hfill & {0.5882} \hfill & { - \,0.9412} \hfill \\ \end{array} } \right]$$
$$\varPi^{w} = \left[ {\begin{array}{*{20}l} { - \,0.7916} \hfill & {0.3332} \hfill & {0.2084} \hfill & {0.25} \hfill \\ {0.2727} \hfill & { - 0.8182} \hfill & {0.1364} \hfill & {0.4091} \hfill \\ {0.1429} \hfill & {0.2857} \hfill & { - \,0.9048} \hfill & {0.4762} \hfill \\ {0.2903} \hfill & {0.1291} \hfill & {0.3545} \hfill & { - \,0.7744} \hfill \\ \end{array} } \right]$$

By substituting the matrix, the corresponding transition probabilities and the positive definite symmetric matrix into the linear matrix inequalities (28), the condition is satisfied, and the multi-energy generation system (27) is stochastically stable, which means the power generation fluctuations of renewable energies in this system can be weakened by adopting the adjustable load control strategy. Figure 1 exhibits the voltage at the connected point of the multi-energy generation system and the grid. Figures 2 and 3 are the output current and output power of wind power. Figures 4 and 5 exhibit the output current and output power of the PV system. Apparently, the current and transmission power vary with the change of solar irradiation and wind speed, which is harmful to the safe and efficient operation of power grid.

Fig. 1
figure 1

Voltage at the connection point

Full size image
Fig. 2
figure 2

Output current of wind power system

Full size image
Fig. 3
figure 3

Output power of wind power system

Full size image
Fig. 4
figure 4

Output current of PV system

Full size image
Fig. 5
figure 5

Output power of PV system

Full size image

According to (15)–(16) and the Markov chain, the probability of solar irradiation and wind speed in each Markov mode can be approximately calculated

$$\begin{aligned} \pi_{\text{s}} & = (0.0301,0.3502,0.4492,0.1705) \\ \pi_{\text{w}} & = (0.2352,0.2229,0.2156,0.3263) \\ \end{aligned}$$

where \(\pi_{\text{s}}\) is the probability of solar irradiation in different Markov modes and \(\pi_{\text{w}}\) is the probability of wind speed in different Markov modes.

According to (17) that the output power expectation of the multi-energy system is 170.1906 kw and the mathematical standard deviation is 38.01 kw, the storage device capacity is set as 50 kWh, which is more than the standard deviation, and kp and ki in formula (25) are 1.5 and 200, respectively. After adopting the adjustable loads control strategy, the current and power transmitted to the grid could remain stable continuously. Figure 6 exhibits the output current of the multi-energy system transmitted to the grid.

Fig. 6
figure 6

Output current of multi-energy generation system with adjustable loads control

Full size image

In [14], a coordination control strategy of battery energy storage is used to smooth the output power fluctuations. In order to make a full and fair comparison, in the same conditions, Fig. 7 shows the output power comparison of multi-energy system by adopting battery energy storage coordination control method and adjustable loads control method when the constraints (24) are satisfied. Obviously, by using the adjustable loads control strategy proposed in this paper, it can get a better performance and more stable results. In addition, the proposed control strategy can maximize the utilization of renewable energies and largely suppress renewable energies power fluctuations.

Fig. 7
figure 7

Output power comparison of multi-energy generation system with different control methods

Full size image

As a matter of fact, the constraints of the flexible loads greatly affect the results of multi-energy generation system. Considering the different constraints scope, satisfying the constraints (24) and not satisfying the power constraint, Fig. 8 shows the comparison result of the output power when the maximum power consumption is reduced to different degrees at 1 s. Figure 9 exhibits the comparison of the output power satisfying the constraints (24) and not satisfying the ramp rate constraint. According to the above results, it can be seen that the maximum power consumption determines whether the excess energy can be consumed, and the regulation ramp rate determines whether the results can quickly converge.

Fig. 8
figure 8

Comparison of the output power satisfying the constraints (24) and not satisfying the power constraint

Full size image
Fig. 9
figure 9

Comparison of the output power satisfying the constraints (24) and not satisfying the ramp rate constraint

Full size image

5 Conclusion

This paper mainly completed the following three tasks: First, the PV generation and wind power generation were analyzed, respectively, and a Markov chain was built according to the cluster results of the solar data and the wind data, and then, a coupled Markov model of the multi-energy generation system with integration of wind and PV was developed. Second, the numerical characteristics of the multi-energy generation system were analyzed and the adjustable control strategy was proposed, and then, the stochastic stability of the system was analyzed by using Markov stability theory. Finally, the rationality and effectiveness of the proposed adjustable loads control strategy were verified by simulation results, and the effect on the regulation results of different constraints scope of the flexible loads was compared.