Design of Intelligent Nonlinear H₂/H_∞ Robust Control Strategy of Diesel Generator Based CPSOGSA Optimization Algorithm

1. Introduction

The diesel generator (DG) is a typical energy and power equipment widely used in the human industrial system [1]. In particular, In areas with weak grids and high cost of long-distance power transmission, and some particular power supply scenarios (such as military emergency power supply and disaster relief power supply), traditional large power grids cannot provide the adequate and timely power supply [2,3,4]. In this case, DGs and their systems can work well.

The system consisting of DGs and the various types of loads connected can be considered a micro-grid (MG) [5]. The MG is an effective way to accelerate renewable energy development and solve the shortcomings of large grids [6,7]. The complementation of photovoltaic panels, wind turbines and DGs in MG is one of the development trends of power generation and distribution in the future power system [8]. With the continuous improvement and perfection of DG technology, more and more people will also benefit from it, especially far away from the large grid. Therefore, the DG research work has a progressive significance for energy supply of human society.

In some MGs, the DG is usually configured as the primary power source to maintain the safe and stable operation of the MG [9,10]. For most MGs, DGs account for a large proportion of their capacity, and their control performance dramatically impacts the MG’s stability [11]. Therefore, the research on the control of the DG is fundamental. However, from the current research status, the comprehensive disturbance characteristics of voltage and frequency of distributed MG are rarely considered in the control research process of conventional DGs [12,13].

The speed control system composed of the diesel engine (DE) and the governor determines the frequency of the synchronous generator, and the frequency stability of the MG system depends on the characteristics of the DE governor system [14,15]. When the load in the MG system changes, the power flow of the system and the output power of each generator will also change suddenly, which will destroy the power balance between the generator and the prime mover. Unbalanced torque occurs, and thus relative motion occurs between the rotors of the generators [16]. The frequency fluctuation of the MG system not only affects users but also deteriorates the operation of the entire system. Therefore, frequency adjustment of the prime mover is necessary [17,18]. Most of the work focuses on optimizing the speed controller. Valenzuela et al. [19] proposed a conventional PI speed controller with a reduced bandwidth to suppress any excitation at the resonant frequency. This slow dynamic performance is also undesirable. Mahmoud et al. [20] analyzed the optimal tuning method of PI/PID controller to improve system performance. The main limitation of this technique is that it is ineffective when the inertia of the driving machine is comparable to that of the loading device [21]. In addition, speed controller optimization using advanced control strategies such as sliding mode controller [22], model predictive control [23], fuzzy logic [24,25], neural networks [26] and optimization algorithm [27]. These proposed solutions show improved dynamic properties and robustness to parameter changes, but these works do not consider the effects of the excitation regulation process. In some extreme cases, these methods may not apply to DG.

We hope that the mechanical power of the prime mover is balanced with the electromagnetic power of the generator. Otherwise, it will cause the rotor of the generator to sway and also cause the electromagnetic transient process and the dynamic process of excitation regulation caused by the voltage and current changes at the generator terminal [28,29]. Since the electrical equipment in the DG system is a mostly inductive load, the load current has a demagnetizing effect on the synchronous generator, and the change of the current size and power factor will cause the change of the terminal voltage of the generator, so the synchronous generator must be equipped with excitation control device [30]. At the same time, in order to improve the reliability of the power supply and generate electricity economically and reasonably, some systems will have several generator sets running in parallel according to different working conditions. In order to make the parallel system run stably, the reactive power between the generators must be reasonably assigned, and this task is also regulated by the excitation control [31,32]. In the study of the excitation control problem, Patel et al. [33] proposed an alternative method that follows the principle of command filtered backstepping to control the terminal voltage and stabilize the DG through the excitation system. In addition, many researchers have proposed to introduce other advanced control techniques into excitation control, such as: decentralized improved cooperative excitation controller [34], excitation controller based on linear matrix inequality (LMI) [35], robust control [36,37], Model Predictive Control [38] for better performance. However, some controllers in these methods are developed based on the linear approximation and bounded operating range of the system, so their performance may degrade significantly when the system is subjected to large disturbances [39].

The dynamic process of DGs in the system includes both electromagnetic dynamic process and dynamic mechanical process, and the two are related to each other [40]. The comprehensive control of speed and excitation can be used as an effective measure to improve the transient stability of the MG system [41]. However, most of the existing DG strategies only consider one of these methods while ignoring the dual control functions of voltage regulation and prime mover speed regulation. In the above studies, stable islanded MGs based on nonlinear control techniques for speed regulation and excitation robust controllers for DGs have been studied [42]. However, the speed regulation and robust excitation controller they designed is not intelligent; that is, the user manually selects the parameters of the robust controller, and no optimization is used to optimize the adjustable parameters.

This paper proposes a robust control method for DG speed and excitation based on the intelligent CPSOGSA, which is used in island MG with multiple load types. The CPSOGSA optimizes the parameters of the robust controller, and the CPSOGSA is introduced for the robust control of DG speed and excitation. Simulation results show that the method effectively suppressed oscillations under severe disturbances and uncertain conditions. The main innovations of this paper can be divided into the following categories:

-

Based on the theory of direct feedback linearization, a nonlinear speed and excitation robust controller of a DG is designed;

-

Applying the intelligent CPSOGSA to optimize the dynamic output function parameters of the robust controller, thus introducing the CPSOGSA for this problem solving multi-objective mixed $H_2/H_{\infty }$ robust controller

-

The proposed method effectively suppresses frequency and voltage oscillations under various load disturbances and uncertainties;

-

Excellent damping efficiency, especially low overshoot, steady-state error, and settling time.

Figure 1. Widely used energy power device —- diesel generator.

Figure 1. Widely used energy power device —- diesel generator.

This paper is organized as follows: In Section 2, the isolated MG structure supplied by DG studied in this paper is introduced, and the dynamic mathematical model of the DG is established in detail, respectively. Then, in Section 3, a speed and excitation robust controller is designed for the MG system’s primary power source-DG. Section 4 presents the CPSOGSA and its application to robust controller optimization problems. Simulation studies are conducted in Section 5 to demonstrate the advantages of the new method. Finally, Section 6 concludes this paper.

2. Modeling of diesel generator and its loads

This section briefly introduces the nonlinear mathematical model of DG, which mainly consists of two parts: speed regulation and excitation system of DG.

2.1. Modeling of speed control part

As shown in Figure 2, the diesel in the DG is a kind of power machine, which can convert the heat energy in the diesel fuel into the mechanical energy required by the generator [43]. Since the main driving output torque $M_1$ of the diesel is a pulsating torque [44], it can be decomposed into a mathematical expression in the form of a Fourier series, which can be described as:

$$M_1=M_p+\sum _{f=1}^{\infty }{M}_f^{0}sin(f{\omega }_{0}t+{\phi }_f)$$

(1)

where $M_1$ is the DE pulsating torque; $M_p$ is the average torque of one cycle; $M_f^0$ is the fth harmonic torque amplitude; ${\omega }_0$ is the harmonic torque angular velocity; ${\phi }_f$ is the fth harmonic torque initial phase angle.

As shown in Figure 3, a flywheel with a significant moment of inertia is installed on the main shaft of the diesel to keep the output torque of the diesel within the allowable range , which can make the flywheel rotation unevenness between $0.00330.0067$ [45]. During operation, the enormous flywheel torque ensures that the output torque of the diesel is uniform. Therefore, the output torque’s unevenness can be ignored in the diesel’s regular operation. In other words, the first constant in Eq. (1) can be retained, and the rest of the harmonic torque composed of amplitude, phase, and frequency can be ignored, and Eq. (2) can be derived

$$M_1=M_p$$

(2)

The driving torque and rotating speed characteristics of DE are smooth curves, approximating several linear segments [46]. Divide the driving torque and rotating speed characteristics of the diesel into m sections; then each section can be written as

$$M_{1j}=k_j{n}_d+b_j$$

(3)

where $n_d$ is the DE rotational speed; the constant value of $k_j$ and $b_j$ is determined by the positions of different segments ($j=1,2,\cdots ,m$).

It can be seen from the adjustment characteristics of the diesel that the relationship between the main driving torque $M_1$ of the DE and the output throttle displacement L of the actuator presents a specific linear relationship [48]. At the same time, the driving torque of the diesel has a certain lag, then the adjustment characteristic of the diesel can be expressed as

$$M_1=\frac{M_1^e}{L_e-L_0}(L(t-T_d)-L_0)$$

(4)

where $M_1^e$ is the torque at rated oil supply; $L_e$ is the rated stroke of the DE actuator; $L_0$ is the no-load stroke of the DE actuator; $T_d$ is the DE torque lag time.

Through the above analysis and derivation, the rotating speed characteristic and adjustment characteristic of the integrated the DE can be obtained as the expression of the driving torque $M_1$:

$$M_1=k_j{n}_d+b_j+\frac{M_1^e}{L_e-L_0}(L(t-T_d)-L_0)$$

(5)

When the DE is running at the rated speed, the segment described in Eq. (5) can be described as

$$M_1=k_1{n}_d+b_1+\frac{M_1^e}{L_e-L_0}(L(t-T_d)-L_0)$$

(6)

Define $d_1=b_1-\left(M_1^e{L}_0\right)/\left(L_e-L_0\right)$, then Eq. (6) can be written as

$$M_1=k_1{n}_d+d_1+\frac{M_1^e}{L_e-L_0}L(t-T_d)$$

(7)

The rotational speed, the number of cylinders, and the number of strokes of the diesel determine the torque lag time $T_d$ of the diesel [47]. Generally, the following equation can be used to estimate the value range

$${\displaystyle \frac{15}{n_d}}<T_d<{\displaystyle \frac{15}{n_d}}+{\displaystyle \frac{60\tau }{n_d{l}_s}}$$

(8)

where $\tau$ is the DE stroke coefficient; $l_s$ is the number of DE cylinders.

The diesel parameters studied in this paper are: the value of rated speed $n_d=1500r/min$, the number of cylinders is 16, and the stroke coefficient $\tau =2$. According to Eq. (8), the value range of $T_d$ can be calculated as $0.01s<T_d<0.015$. The speed transition time of the diesel is 2 seconds, and $T_d$ is very small compared with it and can be ignored. Therefore, Eq. (7) can become

$$M_1=k_1{n}_d+d_1+\frac{M_1^e}{L_e-L_0}L$$

(9)

where L is the throttle actuator displacement.

The relationship between the rotational speed of the diesel $n_d$ and the angular velocity of the crank shaft ${\omega }_g$ is

$$n_d={\displaystyle \frac{60{\omega }_g}{2\pi }}$$

(10)

where ${\omega }_g$ is the angular velocity of the crankshaft shaft.

Eq. (9) can be written as:

$$M_1={\displaystyle \frac{60k_1{\omega }_g}{2\pi }}+d_1+\frac{M_1^e}{L_e-L_0}L$$

(11)

The dynamic motion equation of the main shaft of the DG set can be expressed as

$$J{\dot{\omega }}_g=M_1-M_2-M_d$$

(12)

where J is the DE rotational moment of inertia; ${\omega }_g$ is the DE spindle angular velocity; $M_1$ is the DE output shaft torque; $M_2$ is the DE load torque; and $M_d$ is the generator damping torque.

The damping torque $M_d$ is generated by the damping winding of the synchronous generator, and its characteristic is that it is proportional to the electrical angular speed of the generator rotor, which can be calculated by Eq. (13).

$$M_d=Dp{\omega }_g$$

(13)

where D is the damping factor; p is the number of synchronous generator pole pairs.

Substituting Eq. (11) and Eq. (13) into Eq. (12) get

$$J{\dot{\omega }}_g+Dp{\omega }_g={\displaystyle \frac{60k_1{\omega }_g}{2\pi }}+d_1+\frac{M_1^e}{L_e-L_0}L-M_2$$

(14)

General variables are marked per unit value in the research and analysis of electric power systems. In order to unify the form, Eq. (14) can be standardized as per unit value. This paper converts the mechanical angular speed into electrical angular speed and then standardizes the power, torque, and angular speed to study the speed response of the diesel. This paper defines the rated apparent power $S_B$ of the DG as the reference value of power and defines ${\omega }_0=100\pi$ as the reference value of the angular speed. Therefore, in standardized form, Eq. (14) can be written as

$${\displaystyle \frac{J{\omega }_{g0}^2}{S_B}}\dot{\omega }={\displaystyle \frac{60k_1{\omega }_{g0}^2-2\pi Dp{\omega }_{g0}^2}{2\pi S_B}}\omega +{\displaystyle \frac{d_1{\omega }_{g0}}{S_B}}+{\displaystyle \frac{a{\omega }_{g0}}{S_B}}L-M_2$$

(15)

where ${\omega }_{g0}$ is $100\pi$.

Ignoring the no-load torque and the winding loss, $M_2$ is equal to the output torque of the synchronous generator. In the stability analysis or control system design, for the convenience of calculation, or it is considered that the angular speed $\omega$ does not change much during the transient process, approximately equal to the synchronous speed, that is, $\omega \approx 1$. Therefore, Eq. (15) can be derived as

$${\displaystyle \frac{J{\omega }_{g0}^2}{S_B}}\dot{\omega }={\displaystyle \frac{60k_1{\omega }_{g0}^2-2\pi Dp{\omega }_{g0}^2}{2\pi S_B}}\omega +{\displaystyle \frac{d_1{\omega }_{g0}}{S_B}}+{\displaystyle \frac{a{\omega }_{g0}}{S_B}}L-P_e$$

(16)

The output power $P_e$ of the salient-pole synchronous generator can be described as

$$P_e={\displaystyle \frac{E_q^{\prime }U}{X_d^{\prime }}}sin\delta +{\displaystyle \frac{U^2}{2}}{\displaystyle \frac{X_d^{\prime }-X_q}{X_d^{\prime }{X}_q}}sin\left(2\delta \right)$$

(17)

where $E_q^{\prime }$ is the q-axis transient potential; U is the DG terminal voltage; X is the generator winding reactance.

The relationship between generator rotor angle $\delta$ and electrical angular velocity $\omega$ be described as

$$\dot{\delta }=(\omega -1){\omega }_0$$

(18)

where ${\omega }_{g0}$ is the diesel generator idle angular speed.

From Eq. (16), (17) and (18), the mathematical model of the electrical transient process of DGs is

$$\left\{\begin{array}{cc}\hfill \dot{\delta }=& (\omega -1){\omega }_0\hfill \\ \hfill \dot{\omega }=& \frac{T_{\mathrm{b}}}{T_{\mathrm{a}}}\omega +\frac{1}{T_{\mathrm{a}}}{c}_1+\frac{c_2}{T_{\mathrm{a}}}L-\hfill \\ \hfill & \frac{1}{T_{\mathrm{a}}}\frac{E_{\mathrm{q}}^{\prime }U}{X_{\mathrm{d}}^{\prime }}sin\delta -\frac{1}{T_{\mathrm{a}}}\frac{U^2}{2}\frac{X_{\mathrm{d}}^{\prime }-X_{\mathrm{q}}}{X_{\mathrm{d}}^{\prime }{X}_{\mathrm{q}}}sin\left(2\delta \right)\hfill \end{array}\right.$$

(19)

in the above, the definitions of $T_a$, $T_b$, $c_1$, $c_2$, and a are as follows

$$\left\{\begin{array}{cc}\hfill & T_a={\displaystyle \frac{J{\omega }_{g0}^2}{S_B}}\hfill \\ \hfill & T_b={\displaystyle \frac{60k_1{\omega }_{g0}^2-2\pi Dp{\omega }_{g0}^2}{2\pi S_B}}\hfill \\ \hfill & c_1={\displaystyle \frac{d_1{\omega }_{g0}}{S_B}}\hfill \\ \hfill & c_2={\displaystyle \frac{a{\omega }_{g0}}{S_B}}\hfill \\ \hfill & a=\frac{M_1^e}{L_e-L_0}\hfill \end{array}\right.$$

2.2. Modeling of excitation control part

During the operation of DGs, excitation control is the most basic and indispensable. The synchronous generator is the control object of the excitation control system [48]. Studying the dynamic characteristics of the excitation system is inseparable from the analysis of the dynamic characteristics of the synchronous generator. The mathematical model of the electromagnetic transient process of the DG includes the stator voltage balance equation of the generator and the transient electromagnetic equation of each winding of the rotor [49]. The standard form of the transient electromagnetic equation of the synchronous generator can be described as

$$\left\{\begin{array}{cc}\hfill & {\dot{E}}_q^{\prime }={\displaystyle \frac{1}{T_{d0}}}{E}_{fd}-{\displaystyle \frac{1}{T_{d0}}}{E}_q^{\prime }-{\displaystyle \frac{X_d-X_d^{\prime }}{T_{d0}}}\hfill \\ \hfill & {\dot{E}}_q^{\prime \prime }={\displaystyle \frac{1}{T_{d0}}}{\displaystyle \frac{X_d^{\prime \prime }-X_l}{X_d^{\prime }-X_l}}{E}_{fd}+\left({\displaystyle \frac{1}{T_{d0}^{\prime \prime }}}-{\displaystyle \frac{X_d^{\prime \prime }-X_l}{X_d^{\prime }-X_1}}{\displaystyle \frac{1}{T_{d0}}}\right)E_q^{\prime }\hfill \\ & -{\displaystyle \frac{1}{T_{d0}^{\prime \prime }}}{E}_q^{\prime \prime }-\left({\displaystyle \frac{X_d^{\prime }-X_d^{\prime \prime }}{T_{d0}^{\prime \prime }}}+{\displaystyle \frac{X_d^{\prime \prime }-X_l}{X_d^{\prime }-X_l}}{\displaystyle \frac{X_d-X_d^{\prime }}{T_{d0}}}{I}_d\right)\hfill \\ \hfill & {\dot{E}}_d^{\prime \prime }=-{\displaystyle \frac{1}{T_{q0}^{\prime \prime }}}{E}_d^{\prime \prime }+\frac{X_q-X_q^{\prime \prime }}{T_{q0}^{\prime \prime }}{I}_q\hfill \\ \hfill & U_d=-R{I}_d+\omega X_q^{\prime \prime }{I}_q+\omega E_d^{\prime \prime }\hfill \\ \hfill & U_q=-R{I}_q-\omega X_d^{\prime \prime }{I}_d+\omega E_q^{\prime \prime }\hfill \\ \hfill & U=\sqrt{U_d^2+U_q^2}\hfill \end{array}\right.$$

(20)

where U is the stator winding terminal voltage; $U_d$ and $U_q$ are the d-axis and q-axis components of the stator winding terminal voltage; R is the stator winding resistance; X is the winding reactance; I is the winding current; T is the winding time constant; $E_d^{\prime \prime }$ is the d-axis subtransitory potential; $E_q^{\prime }$ is the q-axis transient potential; $E_q^{\prime \prime }$ is the q-axis subtransitory potential; $E_{df}$ is the excitation winding voltage.

Combining Eq. (19) and Eq. (20), a nonlinear mathematical model of the DG is obtained. This model reflects the interaction and mutual influence of speed and voltage and more accurately describes the change law of speed and voltage.

3. Design of nonliear $H_2/H_{\infty }$ synthetic controller

3.1. Multi-objective state-feedback theory

The control system shown in Figure 4 defines the closed-loop transfer function $T_{\infty }\left(s\right)$ from disturbance w to $H_{\infty }$ control performance output vector $z_{\infty }$, and the closed-loop transfer function $T_2\left(s\right)$ from disturbance w to $H_2$ control performance output vector $z_2$. The control goal is to design a state feedback control law $u=Kx$ to stabilize the closed-loop system progressively. In addition, it needs to meet:

-

Keep the RMS gain ($H_{\infty }$ norm) of $T_{\infty }\left(s\right)$ below a certain specified value ${\gamma }_0>0$;

-

Maintain the $H_2$ norm (LQG cost) of $T_2\left(s\right)$ below a certain specified value ${\upsilon }_0>0$;

-

Minimized form of $H_2/H_{\infty }$ trade-off standard $\alpha T_2\left(s\right)+\beta T_{\infty }\left(s\right)$;

-

Place the closed-loop pole in the designated area of the open left half-plane.

The mixed $H_2/H_{\infty }$ standard considers both the interference suppression aspect (RMS gain from w to $z_{\infty }$) and the LQG aspect ($H_2$ norm from w to $z_2$). In addition, the closed-loop pole can be forced into a specific part of the stable half-plane to obtain a well-damped transient response. Such control problems are called $H_2/H_{\infty }$ control design problems [49,50].

The linear system shown in Figure (Figure 4) can be expressed as:

$$\left\{\begin{array}{cc}\hfill & \dot{x}=\mathit{A}x+{\mathit{B}}_{1}w+{\mathit{B}}_{2}u\hfill \\ \hfill & z_{\infty }={\mathit{C}}_{1}x+{\mathit{D}}_{11}w+{\mathit{D}}_{12}u\hfill \\ \hfill & z_2={\mathit{C}}_{2}x+{\mathit{D}}_{21}w+{\mathit{D}}_{22}u\hfill \end{array}\right.$$

(21)

where the first term in Eq. (21) is a dynamic linear state-space differential equation, and the second and third terms are its dynamic performance evaluation signals.

The closed-loop system formed after adding the controller $\mathit{K}x$ can be described in the form of state-space as:

$$\left\{\begin{array}{cc}\hfill & \dot{x}={\mathit{A}}_{cl}x+{\mathit{B}}_{1}w\hfill \\ \hfill & z_{\infty }={\mathit{C}}_{cl1}x+{\mathit{D}}_{11}w\hfill \\ \hfill & z_2={\mathit{C}}_{cl2}x+{\mathit{D}}_{21}w\hfill \end{array}\right.$$

(22)

where ${\mathit{A}}_{cl}=\mathit{A}+{\mathit{B}}_2\mathit{K}$; ${\mathit{C}}_{cl1}={\mathit{C}}_1+{\mathit{D}}_{12}\mathit{K}$; ${\mathit{C}}_{cl2}={\mathit{C}}_2+{\mathit{D}}_{22}\mathit{K}$.

Taken separately, the two design goals have the following LMI formulation:

$H_{\infty }$ performance: the closed-loop RMS gain from w to $z_{\infty }$ does not exceed $\gamma$ if and only if there exists a symmetric matrix $X_{\infty }$ such that

$$\left\{\begin{array}{cc}\hfill & \left[\begin{array}{ccc}{\mathit{A}}_{cl}{\mathit{X}}_{\infty }+{\mathit{X}}_{\infty }{\mathit{A}}_{cl}^T& {\mathit{B}}_1& {\mathit{X}}_{\infty }{\mathit{C}}_{cl1}^T\\ {\mathit{B}}_1^T& -\gamma \mathit{I}& {\mathit{D}}_{11}^T\\ {\mathit{C}}_{cl1}{\mathit{X}}_{\infty }& {\mathit{D}}_{11}& -\gamma \mathit{I}\end{array}\right]<0\hfill \\ \hfill & {\mathit{X}}_{\infty }>0\hfill \end{array}\right.$$

(23)

$H_2$ performance: the closed-loop $H_2$ norm of $T_2\left(s\right)$ does not exceed $\upsilon$ if there exist two symmetric matrices ${\mathit{X}}_2$ and $\mathit{Q}$ such that

$$\left\{\begin{array}{cc}\hfill & \left[\begin{array}{cc}{\mathit{A}}_{cl}{\mathit{X}}_2+{\mathit{X}}_2{\mathit{A}}_{cl}^T& {\mathit{B}}_1\\ {\mathit{B}}_1^T& -\mathit{I}\end{array}\right]<0\hfill \\ \hfill & \left[\begin{array}{cc}\mathit{Q}& {\mathit{C}}_{cl2}{\mathit{X}}_2\\ {\mathit{X}}_2{\mathit{C}}_{cl2}^T& {\mathit{X}}_2\end{array}\right]>0\hfill \\ \hfill & Trace\left(\mathit{Q}\right)<{\upsilon }^2\hfill \end{array}\right.$$

(24)

The mixed $H_2/H_{\infty }$ performance: The pole of the closed-loop system is located in the left half-open complex plane, and the performance index ${\alpha \parallel T\left(s\right)}_{\infty }{\parallel }_{\infty }^2+\beta {\parallel T_2\left(s\right)\parallel }_2^2$ is minimized.

To solve the mixed objective, the two sets of conditions add up to a non-convex optimization problem with variables $\mathit{Q}$, $\mathit{K}$, ${\mathit{X}}_{\infty }$, and ${\mathit{X}}_2$. In order to facilitate the processing of this problem in the LMI framework, we seek a single Lyapunov matrix $\mathit{X}:={\mathit{X}}_{\infty }={\mathit{X}}_2$ to force the solution of the two goals. With the change of variable $\mathit{Y}:=\mathit{K}\mathit{X}$, this leads to the multi-objective state feedback synthesis problem, which the following sub-optimal LMI formula can express:

Minimize $\alpha {\gamma }^2+\beta Trace\left(\mathit{Q}\right)$ over $\mathit{Y}$, $\mathit{X}$, $\mathit{Q}$, and ${\gamma }^2$ satisfying:

$$\left\{\begin{array}{cc}\hfill & \left[\begin{array}{ccc}\mathit{A}\mathit{X}+\mathit{X}{\mathit{A}}^T+{\mathit{B}}_2\mathit{Y}+{\mathit{Y}}^T{\mathit{B}}_2^T& {\mathit{B}}_1& \mathit{X}{\mathit{C}}_1^T+{\mathit{Y}}^T{\mathit{D}}_{12}^T\\ {\mathit{B}}_1^T& -\gamma \mathit{I}& {\mathit{D}}_{11}^T\\ {\mathit{C}}_1\mathit{X}+{\mathit{D}}_{12}\mathit{Y}& {\mathit{D}}_{11}& -\gamma \mathit{I}\end{array}\right]<0\hfill \\ \hfill & \left[\begin{array}{cc}\mathit{Q}& {\mathit{C}}_2\mathit{X}+{\mathit{D}}_{22}\mathit{Y}\\ \mathit{X}{\mathit{C}}_2^T+{\mathit{Y}}^T{\mathit{D}}_{22}^T& \mathit{X}\end{array}\right]<0\hfill \\ \hfill & Trace\left(\mathit{Q}\right)<{\upsilon }_0^2\hfill \\ \hfill & {\gamma }_0<{\gamma }_0^2\hfill \end{array}\right.$$

(25)

Denoting the optimal solution by (${\mathit{X}}^{*},{\mathit{Y}}^{*},{\mathit{Q}}^{*},{\gamma }^{*}$), the corresponding state-feedback gain is given by:

$${\mathit{K}}^{*}={\mathit{Y}}^{*}{\left({\mathit{X}}^{*}\right)}^{-1}$$

(26)

this gain guarantees the worst-case performances:

$$\left\{\begin{array}{cc}\hfill & {\parallel T\left(s\right)}_{\infty }{\parallel }_{\infty }\le {\gamma }^{*}\hfill \\ \hfill & {\parallel T\left(s\right)}_2{\parallel }_2\le \sqrt{Trace\left({\mathit{Q}}^{*}\right)}\hfill \end{array}\right.$$

(27)

The solution of the mixed $H_2/H_{\infty }$ controller is usually configured with $\alpha =0.5$ and $\beta =0.5$ to obtain ([48]). The robust controller obtained by this method can improve and optimize the control performance of the related system.

3.2. Design of speed controller

The schematic diagram of the DG integrated control system based on the nonlinear $H_2/H_{\infty }$ integrated controller is shown in Figure 5. The DG integrated control system comprises a diesel, generator, nonlinear $H_2/H_{\infty }$ integrated controller, actuator, fuel supply mechanism, and exciter. The nonlinear $H_2/H_{\infty }$ integrated controller consists of two parts, one part is the nonlinear $H_2/H_{\infty }$ speed controller, and the other is the nonlinear $H_2/H_{\infty }$ voltage controller. The nonlinear $H_2/H_{\infty }$ integrated controller should consider the coupling effect of speed and voltage, comprehensively control the speed and voltage, and control the mutual influence between the two to a minimum. Further, improve the stability of the frequency and voltage of the independent MG system.

The differential equation expression of the actuator is:

$$\dot{L}=-{\displaystyle \frac{1}{T_1}}L+{\displaystyle \frac{K_1}{T_1}}{u}_1$$

(28)

where $T_1$ is the actuator time constant; $K_1$ is the actuator gain.

The expression of the differential equation of the exciter is

$${\dot{E}}_{fd}=-{\displaystyle \frac{1}{T_2}}{E}_{fd}+{\displaystyle \frac{K_2}{T_2}}{u}_2$$

(29)

where $T_2$ is the exciter time constant; $K_2$ is the exciter gain.

Combining Eq. (19) and Eq. (28) can get the nonlinear mathematical model of the DG’s speed control system:

$$\left\{\begin{array}{cc}\hfill \dot{\delta }=& (\omega -1){\omega }_0\hfill \\ \hfill \dot{L}=& -\frac{1}{T_1}L+\frac{K_1}{T_1}{u}_1\hfill \\ \hfill \dot{\omega }=& \frac{T_b}{T_a}\omega +\frac{1}{T_a}{c}_1+\frac{c_2}{T_a}L-\frac{1}{T_a}\frac{E_q^{\prime }U}{X_d^{\prime }}sin\delta -\hfill \\ \hfill & \frac{1}{T_a}\frac{U^2}{2}\frac{X_d^{\prime }-X_q}{X_d^{\prime }{X}_q}sin\left(2\delta \right)\hfill \end{array}\right.$$

(30)

Eq. (30) has apparent nonlinear terms. In order to be able to apply the above mixed $H_2/H_{\infty }$ control theory, this paper uses the direct feedback linearization method to linearize the nonlinear mathematical model of the diesel generator speed system described in Eq. (30) to obtain its linear model and then obtain the corresponding nonlinear speed control law.

Three variables $X_1$, $X_2$ and $X_3$ are defined in this paper, which can be expressed as:

$$\left\{\begin{array}{cc}\hfill X_1=& \delta \hfill \\ \hfill X_2=& (\omega -1){\omega }_0\hfill \\ \hfill X_3=& {\displaystyle \frac{T_b}{T_a}}\omega +{\displaystyle \frac{1}{T_a}}{c}_1+{\displaystyle \frac{c_2}{T_a}}L-{\displaystyle \frac{1}{T_a}}{\displaystyle \frac{E_q^{\prime }U}{X_d^{\prime }}}sin\delta \hfill \\ \hfill & \frac{1}{T_a}\frac{U^2}{2}\frac{X_d^{\prime }-X_q}{X_d^{\prime }{X}_q}sin\left(2\delta \right)\hfill \end{array}\right.$$

(31)

According to the definitions of $X_1$, $X_2$ and $X_3$, the nonlinear differential equations described by Eq. (30) can be rewritten as:

$$\left\{\begin{array}{cc}\hfill {\dot{X}}_1=& X_2\hfill \\ \hfill {\dot{X}}_2=& {\omega }_0{X}_3+d_{1}W\hfill \\ \hfill {\dot{X}}_3=& \frac{T_b}{T_a}{X}_3+\frac{c_2{K}_1}{T_a{T}_1}{u}_1-\frac{c_2}{T_a{T}_1}L-\hfill \\ \hfill & \frac{E_q^{\prime }U{\omega }_0}{T_a{X}_d^{\prime }}cos\delta (\omega -1)-\hfill \\ \hfill & \frac{U^2{\omega }_0\left(X_d^{\prime }-X_q\right)}{T_a{X}_d^{\prime }{X}_d}cos\left(2\delta \right)(\omega -1)\hfill \end{array}\right.$$

(32)

where $d_{1}W$ is the disturbance signal set when the mixed $H_2/H_{\infty }$ control method is adopted, including equivalent disturbances caused by disturbance torque and modelling errors.

For speed control, this paper defines virtual control variables V as:

$$\begin{array}{cc}\hfill V=& {\displaystyle \frac{c_2{K}_1}{T_a{T}_1}}{u}_1-{\displaystyle \frac{c_2}{T_a{T}_1}}L-{\displaystyle \frac{E_q^{\prime }U{\omega }_0}{T_a{X}_d^{\prime }}}cos\delta (\omega -1)-\hfill \\ & {\displaystyle \frac{U^2{\omega }_0(X_d^{\prime }-X_q)}{T_a{X}_d^{\prime }{X}_q}}cos\left(2\delta \right)(\omega -1)\hfill \end{array}$$

(33)

By introducing virtual control variables V, Eq. (30) is described by linear state-space differential equations as:

$$\dot{X}={\mathit{A}}_{s}X+{\mathit{B}}_{1s}W+{\mathit{B}}_{2s}V$$

(34)

where X, ${\mathit{A}}_s$, ${\mathit{B}}_{1s}$ and ${\mathit{B}}_{2s}$ are defined as follows:

$$\begin{array}{cc}\hfill & X=\left[\begin{array}{c}{X}_1\\ X_2\\ X_3\end{array}\right];{\mathit{A}}_s=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& {\omega }_0\\ 0& 0& {\displaystyle \frac{T_{\mathrm{b}}}{T_{\mathrm{a}}}}\end{array}\right];\hfill \\ \hfill & {\mathit{B}}_{1s}=\left[\begin{array}{c}0\\ d_1\\ 0\end{array}\right];{\mathit{B}}_{2s}=\left[\begin{array}{c}0\hfill \\ 0\hfill \\ 1\hfill \end{array}\right].\hfill \end{array}$$

Add the dynamic performance evaluation signal described by Eq. (35) to this system.

$$\left\{\begin{array}{c}\hfill Z_{\infty s}={\mathit{C}}_{1s}X+{\mathit{D}}_{11s}W+{\mathit{D}}_{12s}V\\ \hfill Z_{2s}={\mathit{C}}_{2s}X+{\mathit{D}}_{21s}W+{\mathit{D}}_{22s}V\end{array}\right.$$

(35)

where

$$\begin{array}{cc}\hfill & {\mathit{C}}_{1s}=\left[\begin{array}{ccc}{q}_{11}& 0& 0\\ 0& q_{12}& 0\\ 0& 0& q_{13}\end{array}\right];{\mathit{D}}_{12s}=\left[\begin{array}{c}0\\ 0\\ r_1\end{array}\right];\hfill \\ \hfill & {\mathit{C}}_{2s}=\left[\begin{array}{ccc}{q}_{21}& 0& 0\\ 0& q_{22}& 0\\ 0& 0& q_{23}\end{array}\right];{\mathit{D}}_{22s}=\left[\begin{array}{c}0\\ 0\\ r_2\end{array}\right];\hfill \\ \hfill & {\mathit{D}}_{11s}={\mathit{D}}_{21s}=\left[\begin{array}{c}0\hfill \\ 0\hfill \\ 0\hfill \end{array}\right].\hfill \end{array}$$

${\mathit{C}}_{1s}$, ${\mathit{C}}_{2s}$, ${\mathit{D}}_{11s}$, ${\mathit{D}}_{12s}$, ${\mathit{D}}_{21s}$, ${\mathit{D}}_{22s}$ are weighting matrices, $q_{ij}>0$ ($i=1,2$; $j=1,2,3$) and $r_i>0$ ($i=1,2$) are weighting coefficients. A controller that meets the optimal performance requirements can be obtained by selecting different weighting coefficients. The performance requirements include independent power system stability, frequency adjustment accuracy, and low energy consumption of the speed control system.

From Eqs. (34) and (35), the augmented controlled object based on the mixed $H_2/H_{\infty }$ control theory can be obtained as

$${\mathit{P}}_s=\left[\begin{array}{ccc}{\mathit{A}}_s& {\mathit{B}}_{1s}& {\mathit{B}}_{2s}\\ {\mathit{C}}_{1s}& {\mathit{D}}_{11s}& {\mathit{D}}_{12s}\\ {\mathit{C}}_{2s}& {\mathit{D}}_{21s}& {\mathit{D}}_{22s}\end{array}\right]$$

(36)

Bring the augmented controlled object ${\mathit{P}}_s$ into the framework shown in Figure 4 and meet the design requirements of the multi-objective $H_2/H_{\infty }$ state-feedback controller at the same time, and then the $H_2/H_{\infty }$ state feedback controller can be obtained as:

$$\begin{array}{cc}\hfill V=& {\mathit{F}}_{s}X\hfill \\ \hfill =& \left[\begin{array}{ccc}{f}_{1s}& f_{2s}& f_{3s}\end{array}\right]\left[\begin{array}{c}{X}_1\\ X_2\\ X_3\end{array}\right]\hfill \\ \hfill =& f_{1s}{X}_1+f_{2s}{X}_2+f_{3s}{X}_3\hfill \end{array}$$

(37)

From Eq. (33) and Eq. (37), the speed control law $u_1$ can be solved as:

$$\begin{array}{cc}\hfill u_1=& {\displaystyle \frac{1}{K_1}}L+{\displaystyle \frac{T_a{T}_1}{c_2{K}_1}}{\mathit{F}}_{s}X+{\displaystyle \frac{E_q^{\prime }U{\omega }_0{T}_1}{c_2{K}_1{X}_d^{\prime }}}cos\delta (\omega -1)\hfill \\ \hfill & +{\displaystyle \frac{U^2{\omega }_0{T}_1(X_d^{\prime }-X_q)}{c_2{K}_1{X}_d^{\prime }{X}_q}}cos\left(2\delta \right)(\omega -1)\hfill \end{array}$$

(38)

Substituting $X_1$, $X_2$, and $X_3$ into Eq. (38), the DE nonlinear $H_2/H_{\infty }$ speed control law in the usable form shown in Eq. (39) can be obtained.

$$\begin{array}{cc}\hfill u_1=& {\displaystyle \frac{1}{K_1}}L+{\displaystyle \frac{T_a{T}_1}{c_2{K}_1}}{f}_{1s}\delta +\frac{T_a{T}_1}{c_2{K}_1}{f}_{2s}(\omega -1){\omega }_0+\hfill \\ \hfill & {\displaystyle \frac{T_b{T}_1}{c_2{K}_1}}{f}_{3s}\omega -{\displaystyle \frac{T_1}{c_2{K}_1}}{\displaystyle \frac{E_q^{\prime }U}{X_d^{\prime }}}{f}_{3s}sin\delta -\hfill \\ \hfill & {\displaystyle \frac{T_1}{c_2{K}_1}}{\displaystyle \frac{U^2}{2}}\frac{X_d^{\prime }-X_q}{X_d^{\prime }{X}_q}{f}_{3s}sin\left(2\delta \right)+{\displaystyle \frac{T_1}{c_2{K}_1}}{f}_{3s}{c}_1+\hfill \\ \hfill & {\displaystyle \frac{T_1}{K_1}}{f}_{3s}L+{\displaystyle \frac{E_q^{\prime }U{\omega }_0{T}_1}{c_2{K}_1{X}_d^{\prime }}}cos\delta (\omega -1)+\hfill \\ \hfill & {\displaystyle \frac{U^2{\omega }_0{T}_1(X_d^{\prime }-X_q)}{c_2{K}_1{X}_d^{\prime }{X}_q}}cos\left(2\delta \right)(\omega -1)\hfill \end{array}$$

(39)

3.3. Design of excitation control controller

Combining Eq. (29) and Eq. (20) of the first three terms can get the nonlinear mathematical model of the synchronous generator voltage regulation system, which is shown in Eq. (40).

$$\left\{\begin{array}{cc}\hfill & {\dot{E}}_{fd}=-{\displaystyle \frac{1}{T_2}}{E}_{fd}+{\displaystyle \frac{K_2}{T_2}}{u}_2\hfill \\ \hfill & {\dot{E}}_q^{\prime }={\displaystyle \frac{1}{T_{d0}}}{E}_{fd}-{\displaystyle \frac{1}{T_{d0}}}{E}_q^{\prime }-{\displaystyle \frac{X_d-X_d^{\prime }}{T_d}}{I}_d\hfill \\ \hfill & \dot{\delta }=(\omega -1){\omega }_0\hfill \\ \hfill & \dot{\omega }={\displaystyle \frac{T_b}{T_a{\omega }_0}}\omega +{\displaystyle \frac{1}{T_a{\omega }_0}}{c}_1+{\displaystyle \frac{c_2}{T_a{\omega }_0}}L-\hfill \\ \hfill & {\displaystyle \frac{1}{T_a{\omega }_0}}{\displaystyle \frac{E_q^{\prime }U}{X_d^{\prime }}}sin\delta -{\displaystyle \frac{1}{T_a{\omega }_0}}{\displaystyle \frac{U^2}{2}}{\displaystyle \frac{X_d^{\prime }-X_q}{X_d^{\prime }{X}_q}}sin\left(2\delta \right)\hfill \end{array}\right.$$

(40)

Since the damper winding of the synchronous generator has a weak influence on the excitation control, the fourth and fifth terms of Eq. (20) can be ignored when designing the voltage controller. Regarding the electromagnetic power $P_e$ as an external disturbance, Eq. (40) becomes

$$\left\{\begin{array}{cc}\hfill {\dot{E}}_{fd}=& -{\displaystyle \frac{1}{T_2}}{E}_{fd}+{\displaystyle \frac{K_2}{T_2}}{u}_2\hfill \\ \hfill {\dot{E}}_q^{\prime }=& {\displaystyle \frac{1}{T_{d0}}}{E}_{fd}-{\displaystyle \frac{1}{T_{d0}}}{E}_q^{\prime }-{\displaystyle \frac{X_d-X_d^{\prime }}{T_{d0}}}{I}_d\hfill \\ \hfill \dot{\delta }=& \omega -1\hfill \\ \hfill \dot{\omega }=& {\displaystyle \frac{T_b}{T_a{\omega }_0}}\omega +{\displaystyle \frac{1}{T_a{\omega }_0}}{c}_1+{\displaystyle \frac{c_2}{T_a{\omega }_0}}L-{\displaystyle \frac{1}{T_a{\omega }_0}}{P}_e\hfill \end{array}\right.$$

(41)

Since $E_q^{\prime }$ is not easy to measure, it is not conducive to designing a state feedback controller. Therefore, the voltage error $\Delta U$ is selected as the state variable, so this paper converts $E{q}^{\prime }$ in Eq. (41) into $\Delta U$ by a constant $c_3$. According to the relationship between variable data:

$$U=E_q^{\prime }-c_3\delta$$

(42)

Through the relationship between U and the voltage $U_0$ at the stator winding end of the generator, it can be derived:

$$E_q^{\prime }=U_0+\Delta U+c_3\delta$$

(43)

According to the relationship described in Eq. (43), the state-space differential equation expressed by Eq. (41) can be rewritten as

$$\left\{\begin{array}{cc}\hfill {\dot{E}}_{fd}=& -{\displaystyle \frac{1}{T_2}}{E}_{fd}+{\displaystyle \frac{K_2}{T_2}}{u}_2\hfill \\ \hfill \Delta \dot{U}=& {\displaystyle \frac{1}{T_{d0}}}{E}_{fd}-{\displaystyle \frac{1}{T_{d0}}}\Delta U-{\displaystyle \frac{1}{T_{d0}}}{c}_3\delta -\hfill \\ \hfill & c_3\omega +(c_3-{\displaystyle \frac{1}{T_{d0}}})-{\displaystyle \frac{X_d-X_d^{\prime }}{T_{d0}}}{I}_d\hfill \\ \hfill \dot{\delta }=& \omega -1\hfill \\ \hfill \dot{\omega }=& {\displaystyle \frac{T_b}{T_a{\omega }_0}}\omega +{\displaystyle \frac{1}{T_a{\omega }_0}}{c}_1+{\displaystyle \frac{c_2}{T_a{\omega }_0}}L-{\displaystyle \frac{1}{T_a{\omega }_0}}{P}_e\hfill \end{array}\right.$$

(44)

Eq. (44) is written as a linear state-space differential equation in the form

$${\dot{X}}^{\prime }={\mathit{A}}_e{X}^{\prime }+{\mathit{B}}_{1e}{W}^{\prime }+{\mathit{B}}_{2e}{u}_2$$

(45)

where

$$\begin{array}{cc}\hfill & X^{\prime }\left(t\right)=\left[\begin{array}{c}{E}_{fd}\\ \Delta U\\ \delta \\ \omega \end{array}\right];W^{\prime }=\left[\begin{array}{c}{I}_d\\ L\\ P_e\\ 1\end{array}\right];{\mathit{B}}_{2e}=\left[\begin{array}{c}{\displaystyle \frac{K_2}{T_2}}\\ 0\\ 0\\ 0\end{array}\right];\hfill \\ \hfill & {\mathit{A}}_e=\left[\begin{array}{cccc}-{\displaystyle \frac{1}{T_2}}& 0& 0& 0\\ {\displaystyle \frac{1}{T_{\mathrm{d}0}}}& -{\displaystyle \frac{1}{T_{\mathrm{d}0}}}& -{\displaystyle \frac{c_3}{T_{\mathrm{d}0}}}& -c_3\\ 0& 0& 0& 1\\ 0& 0& 0& {\displaystyle \frac{T_{\mathrm{b}}}{T_{\mathrm{a}}{\omega }_0}}\end{array}\right];\hfill \\ \hfill & {\mathit{B}}_{1e}=\left[\begin{array}{cccc}0& 0& 0& 0\\ X_{\mathrm{d}}-X_{\mathrm{d}}^{\prime }& 0& 0& c_3-{\displaystyle \frac{1}{T_{\mathrm{d}0}}}\\ -{\displaystyle \frac{T_{\mathrm{d}0}}{0}}& 0& 0& -1\\ 0& {\displaystyle \frac{c_2}{T_{\mathrm{a}}{\omega }_0}}& -{\displaystyle \frac{1}{T_{\mathrm{a}}{\omega }_0}}& {\displaystyle \frac{c_1}{T_{\mathrm{a}}{\omega }_0}}\end{array}\right].\hfill \end{array}$$

Similarly, the dynamic performance evaluation signal of the system described in Eq. (45) is defined as

$$\left\{\begin{array}{c}{Z}_{\infty e}={\mathit{C}}_{1e}{X}^{\prime }+{\mathit{D}}_{11e}{W}^{\prime }+{\mathit{D}}_{12e}{u}_2\hfill \\ Z_{2e}={\mathit{C}}_{2e}{X}^{\prime }+{\mathit{D}}_{21e}{W}^{\prime }+{\mathit{D}}_{22e}{u}_2\hfill \end{array}\right.$$

(46)

where

$$\begin{array}{cc}\hfill & {\mathit{C}}_{1e}=\left[\begin{array}{cccc}{q}_{14}& 0& 0& 0\\ 0& q_{15}& 0& 0\\ 0& 0& q_{16}& 0\\ 0& 0& 0& q_{17}\end{array}\right];{\mathit{D}}_{12e}=\left[\begin{array}{c}{r}_3\hfill \\ r_4\hfill \\ r_5\hfill \\ r_6\hfill \end{array}\right];\hfill \\ \hfill & {\mathit{C}}_{2e}=\left[\begin{array}{cccc}{q}_{24}& 0& 0& 0\\ 0& q_{25}& 0& 0\\ 0& 0& q_{26}& 0\\ 0& 0& 0& q_{27}\end{array}\right];{\mathit{D}}_{22e}=\left[\begin{array}{c}{r}_7\\ r_8\\ r_9\\ r_{10}\end{array}\right];\hfill \\ \hfill & {\mathit{D}}_{11e}={\mathit{D}}_{21e}=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right].\hfill \end{array}$$

${\mathit{C}}_{1e}$, ${\mathit{C}}_{2e}$, ${\mathit{D}}_{11e}$, ${\mathit{D}}_{12e}$, ${\mathit{D}}_{21e}$, and ${\mathit{D}}_{22e}$ are weighting matrices, $q_{ij}>0$ ($i=1,2$; $j=4,5,6,7$) and $r_i>0$($i=3,4,\cdots ,10$) are weighting coefficients. The voltage controller that meets the optimal performance is solved by combining different weighting coefficients. The performance includes system stability, voltage regulation accuracy and low energy loss of the excitation system.

From Eqs. (45) and (46), the augmented controlled object ${\mathit{P}}_e$ used to solve the mixed $H_2/H_{\infty }$ robust controller can be obtained, which can be expressed as:

$${\mathit{P}}_e=\left[\begin{array}{ccc}{\mathit{A}}_e& {\mathit{B}}_{1e}& {\mathit{B}}_{2e}\\ {\mathit{C}}_{1e}& {\mathit{D}}_{11e}& {\mathit{D}}_{12e}\\ {\mathit{C}}_{2e}& {\mathit{D}}_{21e}& {\mathit{D}}_{22e}\end{array}\right]$$

(47)

Similarly, the augmented controlled object ${\mathit{P}}_e$ is brought into the solution framework shown in Figure 4, and the voltage $H_2/H_{\infty }$ state feedback controller of the diesel generator is obtained to satisfy the LMI. The nonlinear $H_2/H_{\infty }$ voltage control law of the DG is

$$\begin{array}{cc}\hfill u_2=& {\mathit{F}}_e{X}^{\prime }\hfill \\ \hfill =& \left[\begin{array}{cccc}{f}_{1e}& f_{2e}& f_{3e}& f_{4e}\end{array}\right]\left[\begin{array}{c}{E}_{fd}\\ \Delta U\\ \delta \\ \omega \end{array}\right]\hfill \\ \hfill =& f_{1e}{E}_{fd}+f_{2e}\Delta U+f_{3e}\delta +f_{4e}\omega \hfill \end{array}$$

(48)

Combining the Eq. (39) and Eq. (48) derived above, the diesel generator’s nonlinear robust $H_2/H_{\infty }$ integrated controller is obtained.

$$\left\{\begin{array}{cc}\hfill u_1=& {\displaystyle \frac{1}{K_1}}L+{\displaystyle \frac{T_a{T}_1}{c_2{K}_1}}{f}_{1s}\delta +\frac{T_a{T}_1}{c_2{K}_1}{f}_{2s}(\omega -1){\omega }_0+\hfill \\ \hfill & {\displaystyle \frac{T_b{T}_1}{c_2{K}_1}}{f}_{3s}\omega -{\displaystyle \frac{T_1}{c_2{K}_1}}{\displaystyle \frac{E_q^{\prime }U}{X_d^{\prime }}}{f}_{3s}sin\delta -\hfill \\ \hfill & {\displaystyle \frac{T_1}{c_2{K}_1}}{\displaystyle \frac{U^2}{2}}\frac{X_d^{\prime }-X_q}{X_d^{\prime }{X}_q}{f}_{3s}sin\left(2\delta \right)+{\displaystyle \frac{T_1}{c_2{K}_1}}{f}_{3s}{c}_1+\hfill \\ \hfill & {\displaystyle \frac{T_1}{K_1}}{f}_{3s}L+{\displaystyle \frac{E_q^{\prime }U{\omega }_0{T}_1}{c_2{K}_1{X}_d^{\prime }}}cos\delta (\omega -1)+\hfill \\ \hfill & {\displaystyle \frac{U^2{\omega }_0{T}_1(X_d^{\prime }-X_q)}{c_2{K}_1{X}_d^{\prime }{X}_q}}cos\left(2\delta \right)(\omega -1)\hfill \\ \hfill u_2=& f_{1e}{E}_{fd}+f_{2e}\Delta U+f_{3e}\delta +f_{4e}\omega \hfill \end{array}\right.$$

(49)

It can be found from Eq. (49) that the designed nonlinear control law considers the coupling effect of the DG’s speed and voltage and comprehensively controls the two, which can improve the dynamic accuracy of the DG’s speed and voltage. Further, improve the stability of the MG system connected to the DG.

4. Artificial hybrid PSOGSA with chaotic maps approach and its application to multi-objective $H_2/H_{\infty }$ robust problem

4.1. The proposed CPSOGSA algorithm

The artificial intelligence CPSOGSA algorithm is a hybrid particle swarm optimization and gravitational search algorithm combined with chaotic mapping [51]. Besides, the CPSOGSA is superior and more robust than other methods (such as MSA, GWO, WOA and so on) [52]. Therefore, this paper also adopts CPSOGSA to optimize the robust controller. This subsection provides the necessary description of the algorithm.

4.1.1. Particle swarm optimization

The PSO consists of three crucial operators, which are the inertia weight vector $w_p$, the personal particle best (pbest), and the global particle best (gbest). The inertia weight gives diversification (exploration) capability, and <pbest, gbest> provides intensification (exploitation) power to PSO, respectively [53]. The velocity $v_i\left(k\right)$ and position $p_i\left(k\right)$ of each particle in the PSO algorithm are updated as follows

$$\left\{\begin{array}{cc}\hfill & v_i(k+1)=w_p{v}_i\left(k\right)+c_1\times \mathrm{rand}\times \left({\mathrm{pbest}}_i-p_i\left(k\right)\right)\hfill \\ & +c_2\times \mathrm{rand}\times \left({\mathrm{gbest}}_i-p_i\left(k\right)\right)\hfill \\ \hfill & p_i(k+1)=p_i\left(k\right)+v_i(k+1)\hfill \end{array}\right.$$

(50)

where $c_1$ and $c_2$ are the acceleration coefficient; rand is a uniform random variable between O and 1.

4.1.2. Gravitational search algorithm

GSA is an optimization algorithm inspired by Newton’s theories of gravity and motion. The operation steps of GSA are briefly summarized as follows [54].

Positions of N agents constituting the initialization population $P_i$ in GSA are listed in Eq. (51). The best and worst fitness values are computed for each iteration. Best and worst fitness values are defined as in Eq. (52) for the problem.

$$P_i=\left(p_1^1,\dots ,p_i^d,\dots ,p_i^n\right),i=1,2,\dots ,N$$

(51)

where $p_i^d$ is the ith agent of the dth dimensions.

$$\begin{array}{cc}\hfill & best\left(k\right)=min{\mathrm{fit}}_i\left(k\right)i\in \{1,2,\dots ,N\}\hfill \\ \hfill & worst\left(k\right)=max{\mathrm{fit}}_i\left(k\right)i\in \{1,2,\dots ,N\}\hfill \end{array}$$

(52)

where ${\mathrm{fit}}_i\left(k\right)$ denotes the fitness value of the ith individual.

The gravitational constant $G\left(k\right)$ for each iteration and inertial masses are calculated as follows

$$\left\{\begin{array}{cc}\hfill & G\left(k\right)=G_0{e}^{\left(-a\frac{k}{K}\right)}\hfill \\ \hfill & M_{ai}=M_{pi}=M_{ii}=M_{i}i=1,2,3,\dots ,N\hfill \\ \hfill & m_i\left(k\right)=\frac{{\mathrm{fit}}_i\left(k\right)-\mathrm{worst}\left(k\right)}{\mathrm{best}\left(k\right)-\mathrm{worst}\left(k\right)}\hfill \\ \hfill & M_i\left(k\right)=\frac{m_i\left(k\right)}{{\sum }_{j=1}^N{m}_j\left(k\right)}\hfill \end{array}\right.$$

(53)

where $G_0$ is the initial value of gravitational constant; k and K are the index and the total number of iterations; $M_{ai}$ is the active mass at kth iteration for the ith agent; $M_{pi}$ is the passive mass at kth iteration for the ith agent; $M_{ii}$ is the the inertial mass at kth iteration for the ith agent; $M_i$ is the mass at kth iteration for the ith agent.

The force $F_{ij}^d\left(k\right)$ acting on the ith mass from the jth mass, the total force $F_i^d\left(k\right)$, and the acceleration value ${\sigma }_i^d\left(k\right)$ at iteration k are calculated as follows

$$\left\{\begin{array}{cc}\hfill & F_{ij}^d\left(k\right)=G\left(k\right)\frac{M_{pi}\left(k\right)\times M_{aj}\left(k\right)}{R_{ij}\left(k\right)+\epsilon }\left(p_j^d\left(k\right)-p_i^d\left(k\right)\right)\hfill \\ \hfill & F_i^d\left(k\right)=\sum _{j\in \left(k\right)\mathrm{best},j\ne i}{\mathrm{rand}}_j{F}_{ij}^d\left(k\right)\hfill \\ \hfill & {\sigma }_i^d\left(k\right)=\frac{F_i^d\left(k\right)}{M_{ii}\left(k\right)}\hfill \end{array}\right.$$

(54)

where $R_{ij}$ is the euclidean distance between the ith mass and the jth mass.

The velocity $V_i^d\left(k\right)$ and position $p_i^d\left(k\right)$ of the ith agent at the dth dimension of the next iteration are updated as follows

$$\left\{\begin{array}{cc}\hfill & V_i^d(k+1)=\mathrm{rand}\times V_i^d\left(k\right)+{\sigma }_i^d\left(k\right)\hfill \\ \hfill & p_i^d(k+1)=p_i^d\left(k\right)+V_i^d(k+1)\hfill \end{array}\right.$$

(55)

4.1.3. Modified hybrid particle swarm optimization and gravitational search algorithm with chaotic maps

The hybrid PSOGSA combines the advantages of PSO and GSA algorithms. The algorithm can be explained a follows

$$\left\{\begin{array}{cc}\hfill & V_i(k+1)=w_p\times V_i\left(k\right)+c_1^{\prime }\times \mathrm{rand}\times {\sigma }_i\left(k\right)\hfill \\ & +c_2^{\prime }\times \mathrm{rand}\times \left(\mathrm{gbest}-p_i\left(k\right)\right)\hfill \\ \hfill & p_i(k+1)=p_i\left(k\right)+V_i(k+1)\hfill \end{array}\right.$$

(56)

When the chaotic map is used in the optimization algorithm, the optimization algorithm can quickly converge to the optimal solution and get rid of the optimal local solution [55]. Therefore, the one-dimensional and irreversible chaotic mapping method with this advantage is applied to the PSOGSA optimization algorithm so that the performance of PSOGSA can be further improved. The chaotic map equation gives these characteristics of the chaotic map. In [56], ten chaotic mapping methods are given to improve the performance of the hybrid PSOGSA algorithm, and the change curve of the chaotic map with the number of iterations is shown in Figure 6. From the figure, it can be found that there is no random value in any chaotic map.

In the proposed CPSOGSA algorithm, the chaotic map method is used to improve the local search ability of GSA. Eqs. (54) calculates the total force, while a random number between 0 and 1 is used as a weight. In other words, this parameter will impact the local search ability of the algorithm. If this situation is considered, using chaotic map methods instead of the random number could increase the convergence ability of the PSOGSA algorithm in an optimization process. Although each chaotic map follows the order in Table, Eq. (57) is applied instead of using a random parameter for calculating the total force in Eqs. (54). The application of chaotic map methods to the total force value is shown in Eq. (57). After the total force value is calculated, the acceleration value is determined using Eqs. (53), while Eqs. (56) updates the velocities and positions of all agents in the population.

$$F_i^d\left(k\right)=\sum _{j\in k\mathrm{best},j\ne i}C\left(k\right)F_{ij}^d\left(k\right)$$

(57)

where $C\left(k\right)$ is the chaotic sequence constant for the kth iteration.

4.2. Applying CPSOGSA to multi-objective robust $H_2/H_{\infty }$ problem

In subsection 3.2, the dynamic performance evaluation function and its weighting matrix for solving the robust controller of DG speed are proposed. The value of the weighting coefficient in these weighting matrices dramatically influences the performance of the speed controller. Therefore, These parameters need to be sufficient, and there are ten weighting coefficients for solving the speed robust controller, including the values of $q_{ij}>0(i=1,2;j=1,2,3)$, $r_i(i=1,2)$. This paper uses the CPSOGSA algorithm to optimize and adjust these parameters. In optimizing the speed controller, the objective optimization function of CPSOGSA is $O{F}_1$. The function described by Eq. (58) is used to calculate the fitness value passed in the optimization process, aiming to minimize the absolute value of the change in the absolute value of the DG speed during this period.

$$O{F}_1=\int {\left|\Delta \omega \right|}^{2}dt$$

(58)

where $\Delta \omega$ is the speed deviation of DG.

Similarly, CPSOGSA is used to adjust the dynamic performance evaluation function of the DG excitation controller proposed in subsection 3.3 and the weighting coefficients in the weighting matrix. For the solution of the excitation controller, the weighting coefficients include $q_{ij}(i=1,2;j=4,5,6,7)$, $r_i(i=3,4,\cdots ,10)$. In the optimization process, the objective optimization function of CPSOGSA is $O{F}_2$. The function described by Eq. (59) is used to calculate the fitness value of the excitation controller during the optimization process, aiming to minimize the absolute value of the DG voltage in time.

$$O{F}_2=\int {\left|\Delta U\right|}^{2}dt$$

(59)

where $\Delta U$ is the voltage deviation of DG.

In conclusion, these parameters need to be adequately tuned. The number of parameters to be adjusted is 16. These parameters are optimally adjusted using the CPSOGSA algorithm, whose overall optimal objective function is considered to be Eq. (60) as the weighted sum of $O{F}_1$ and $O{F}_2$ given by Eqs. (58) and (58). Therefore, the optimization objective can be described as

$$OF=min\left\{O{F}_1+O{F}_2\right\}$$

(60)

5. Numerical study

5.1. Simulation parameters

The MG system composed of DG is shown in Figure 7. Among them, the rated power of the DE is 1250kW, the rated speed is 1500r/min, the rated torque is 11.9kNm, and the moment of inertia is 71.82; while the rated voltage of the generator is 380V, the rated current is 2310A, the power factor is 0.8, and the rated frequency is 50Hz; in addition, there are static and dynamic loads in this MG. The static load is mainly composed of electric heaters and lighting equipment, which accounts for a small proportion of the total loads; the dynamic loads refers to PMSM and driving devices, such as washing machines, water pumps, and fans, which account for a large proportion of the total load. Some parameters are also used for simulation, as shown in Table 1.

The selection of the relevant coefficients in the weighting function is a difficult point in multi-objective $H_2/H_{\infty }$ control and must be selected through repeated trials. After each trial selection, use the LMI toolbox to find the state feedback coefficient, substitute it into the simulation model, and then conduct the characteristic test to obtain the best comprehensive performance index. Simulations are carried out in MATLAB/Simulink and its programming environment interface. The parameters generated by the CPSOGSA algorithm are brought into the weighting function, and the state feedback control matrix ${\mathit{F}}_s$ of the speed regulation part of the DG and the state feedback control matrix ${\mathit{F}}_e$ of the excitation control part are calculated by the LMI and sent to Simulink. Simulink then iteratively evaluates the performance of the DG based on the adjusted parameters. In general, the flowchart of Figure 8 can be thought of as an optimization process. The number of iterative optimizations in this paper is 100 times.

As shown in Figure 9, under the optimization objective described in Eq. (60), the fitness value (OF value) change curve in different algorithm iterations for 100 times. The optimization range of each parameter is $[11e-3]$. It is not difficult to find from the figure that the CPSOGSA algorithm only needs about ten iterations to complete the optimization. On the contrary, the traditional algorithms are insufficient in convergence and optimization efficiency. For example, the particle swarm algorithm PSO needs nearly 80 generations to converge, while the gravitational search algorithm falls into the local optimum at the beginning, and after introducing the chaotic sequence, GSA can avoid falling into the local optimum.

For the design of the speed control part of the nonlinear integrated controller based on the mixed $H_2/H_{\infty }$ control theory, after the optimization of the CPSOGSA optimization algorithm, the weighting matrixs in Eq. (35) are

$$\begin{array}{cc}\hfill & {\mathit{C}}_{1s}={\mathit{C}}_{2s}\left[\begin{array}{ccc}0.4788& 0& 0\\ 0& 0.0015& 0\\ 0& 0& 0.1105\end{array}\right]\hfill \\ \hfill & {\mathit{D}}_{12s}={\mathit{D}}_{22s}\left[\begin{array}{c}0\\ 0\\ 0.0011\end{array}\right]\hfill \end{array}$$

Combining these matrices and using the Linear Matrix Inequality Toolbox in MATLAB, find the feedback coefficients in Eq. (37) as $f_{1s}=-1.9188$, $f_{2s}=-1.8738$, and $f_{3s}=35.604$.

Similarly, for the weighting matrix in Eq. (46), These matrix are

$$\begin{array}{cc}\hfill & {\mathit{C}}_{1e}={\mathit{C}}_{2e}=\left[\begin{array}{cccc}0.7817& 0& 0& 0\\ 0& 3.824\times {10}^3& 0& 0\\ 0& 0& 0.0256& 0\\ 0& 0& 0& 0.1270\end{array}\right]\hfill \\ \hfill & {\mathit{D}}_{12e}={\mathit{D}}_{22e}=\left[\begin{array}{c}0.1038\hfill \\ 0.1143\hfill \\ 0.1043\hfill \\ 0.1417\hfill \end{array}\right]\hfill \end{array}$$

Using the LMI toolbox in Matlab, get the coefficients of state-feedback controller in Eq. (48): $f_{1e}=-2.227$, $f_{2e}=-7.516\times {10}^3$, $f_{3e}-8.328\times {10}^{-4}$ and $f_{4e}=-9.495\times {10}^{-5}$.

5.2. Simulation results

In MATLAB/SIMULINK, the simulation experiment of a diesel generator with the mixed $H_2/H_{\infty }$ robust integrated controller is established to complete the effect verification of the controller described by Eq. (49). The relevant simulation parameters are shown in Table 1. In the comparative simulation process, the diesel generator adopts the optimized PID parameters $K_p=20$, $K_i=0.02$, $K_d=5$.

Firstly, the comparison simulation experiment of diesel generator no-load starting characteristics is carried out. When no load is connected, the diesel generator start-up process is controlled by the controller and PID controller proposed in this paper, respectively. The simulation results are shown in Figure 10 and Figure 11. It can be seen from the speed curve diagram in Figure 10: from 0 to the establishment of the speed of 1 p.u, the proposed controller is about 3s faster than the PID controller. From the voltage change curve in Figure 11, it can be seen that from 0 to 1 p.u of the machine terminal voltage, our controller is 10s faster than the PID controller and has better start-up performance.

Then the computer simulation of the designed system with sudden static loads is carried out. Figure 12 and Figure 13 show the dynamic characteristic curves of speed and voltage of the system after sudden static loads, and the two variables are per unit value. Among them, the blue curve represents the action result of the intelligent nonlinear $H_2/H_{\infty }$ integrated controller, and the red curve represents the action result of the conventional PID controller.

From the simulation results, the dynamic performance of the proposed controller is excellent. The dynamic regulation rate is minimal, the stabilization time is about 2 seconds, the dynamic voltage regulation rate is only about 0.05, and the stabilization time is 4 seconds. Then, when the conventional PID controller is applied, the dynamic regulation rate of the system is significant, and the stabilization time takes 5 seconds. The dynamic regulation rate of the system exceeds 0.15, and the stabilization time takes a long time. It can be seen through computer simulation that after using the proposed controller, the rate and voltage rate are both decreased, and the stabilization time is shortened.

Finally, a simulation is performed that introduces a simulation of a dynamic loads such as PMSM in the system. The dynamic loads greatly influences the grid voltage, especially when the grid is unloaded. The simulation runs for 55 seconds to connect to a PMSM with a size of 90 kW. The simulation results are shown in Figure 14 and Figure 15. From the situation of the speed change in Figure 14, when the controller proposed in this paper is used, the system shows very little oscillation, while when the PID controller is used, the number of speed oscillations of the system exceeds five times. Figure 15 shows the response process of the terminal voltage under this working condition, and the controller proposed in this paper also shows better performance.

In this paper, the intelligent nonlinear integrated controller effectively improves the dynamic accuracy of speed control and excitation control. The improvement of dynamic accuracy weakens the mutual influence and interaction between speed and voltage and solves the comprehensive control problem of speed and voltage. The method improves the stability of the MG system.

6. Conclusion

This paper first establishes a DG’s nonlinear dynamic model. Then, the DG’s speed and excitation controller are designed, respectively. The direct feedback linearization method and the multi-objective robust control theory are applied in the design process. Given the difficulty in determining the weighting parameters in the design process, this paper introduces the artificial intelligence CPSOGSA optimization algorithm to optimize the relevant weighting parameters that affect the controller’s performance. In order to verify the effectiveness of the proposed controller in improving the MG system supplied by DG, a MG model driven by a DG was built in Simulink. The relevant simulation results show that, compared with the traditional PID controller, the proposed controller can effectively improve the dynamic accuracy of the system and the ability to suppress disturbances and reduce the fluctuation of the MG system in the process of sudden load change and operating condition switching. The stability of the independent MG is improved.