*Article* **Current Control of the Permanent-Magnet Synchronous Generator Using Interval Type-2 T-S Fuzzy Systems**

#### **Yuan-Chih Chang \*, Chi-Ting Tsai and Yong-Lin Lu**

Department of Electrical Engineering and Advanced Institute of Manufacturing with High-tech Innovations, National Chung Cheng University, Chiayi 62102, Taiwan

**\*** Correspondence: ycchang@ccu.edu.tw; Tel.: +886-5-272-9108

Received: 10 July 2019; Accepted: 29 July 2019; Published: 31 July 2019

**Abstract:** The current control of the permanent-magnet synchronous generator (PMSG) using an interval type-2 (IT2) Takagi-Sugeno (T-S) fuzzy systems is designed and implemented. PMSG is an energy conversion unit widely used in wind energy generation systems and energy storage systems. Its performance is determined by the current control approach. IT2 T-S fuzzy systems are implemented to deal with the nonlinearity of a PMSG system in this paper. First, the IT2 T-S fuzzy model of a PMSG is obtained. Second, the IT2 T-S fuzzy controller is designed based on the concept of parallel distributed compensation (PDC). Next, the stability analysis can be conducted through the Lyapunov theorem. Accordingly, the stability conditions of the closed-loop system are expressed in Linear Matrix Inequality (LMI) form. The AC power from a PMSG is converted to DC power via a three-phase six-switch full bridge converter. The six-switch full bridge converter is controlled by the proposed IT2 T-S fuzzy controller. The analog-to-digital (ADC) conversion, rotor position calculation and duty ratio determination are digitally accomplished by the microcontroller. Finally, simulation and experimental results verify the performance of the proposed current control.

**Keywords:** permanent-magnet synchronous generator; T-S fuzzy; current control; interval type-2

#### **1. Introduction**

A permanent-magnet synchronous generator (PMSG) [1] is an essential unit implemented to convert energy. Its power density and power conversion efficiency is high and its maintenance cost is low. Consequently, PMSG can be implemented in various utilizations like electric vehicles (EV) and hybrid electric vehicles (HEV) [2,3], home appliances [4], wind generation systems [5], flywheel energy storage systems (FESS) [6] and ultrahigh-speed elevators [7]. The embedded type of the permanent magnet will affect the characteristics of a PMSG. The interior-type [8–10] has high inductance saliency and generates higher torque. However, its torque ripple is also higher. On the contrary, the torque ripple of the surface-mounted type [10,11] is lower and its reluctance torque is nearly absent [12]. In motor design concept, magnet design methods [13,14] are presented to optimize the performance and minimize the torque ripple.

If the d-axis current is controlled at zero, then the electromagnetic torque of a PMSG is proportional to the q-axis current. This characteristic enhances the importance of current controls [15–20] for the PMSG. Since the power generating performance of a PMSG is affected by its winding current, the torque ripple can be eliminated by reducing winding current harmonics [16]. Moreover, the generating capability of a PMSG can be improved via adopting proper and good current control algorithms. The traditionally implemented current control schemes include fixed-frequency control [17,18], hysteresis control [19] and predictive control [20]. The state space equation of a PMSG is nonlinear. Therefore, the Takagi-Sugeno fuzzy [21,22] (T-S fuzzy) system is implemented to design a speed

controller for the permanent-magnet synchronous motor (PMSM) and a current controller for the PMSG. System uncertainty [23] generally exists in nonlinear systems. Therefore, the interval type-2 (IT2) fuzzy logic system [24,25] is proposed to deal with this problem. In this paper, IT2 T-S fuzzy models are implemented in the design of a current controller for the PMSG.

In this study, system configuration of the PMSG is first introduced. Then the dynamic model of a PMSG is conducted. Next, IT2 T-S fuzzy models of a PMSG are established for the design of an IT2 T-S fuzzy current controller. The stability of the IT2 T-S fuzzy control system for PMSG is analyzed using the Lyapunov theorem. The stability conditions of the proposed current controller are expressed in Linear Matrix Inequality (LMI) form. Experimental results, including constant current command tracking, variable current command tracking and computation time of the microcontroller, are demonstrated to verify the performance of the designed IT2 T-S fuzzy current controller.

#### **2. System Configuration and Dynamic Model**

Figure 1 introduces the system configuration of the PMSG based on the IT2 T-S fuzzy control systems. The prime mover of the PMSG consists of a PMSM with a PMSM driver. The input torque of PMSG is provided by the prime mover. The encoder signals on the shaft are utilized to calculate rotor speed and rotor position. Three-phase winding currents are sensed through the analog-to-digital converter (ADC). Then the abc-frame currents are transformed to dq-frame currents. The current commands and dq-frame currents are implemented in the IT2 T-S fuzzy current controller to calculate the dq-frame voltage commands. The corresponding duty ratio of six switches is determined via space-vector pulse-width modulation (SVPWM). All the control and transformation schemes are digitally realized using the microcontroller Renesas RX62T (Renesas Electronics Corp., Tokyo, Japan).

**Figure 1.** System configuration of the PMSG based on theIT2 T-S fuzzy control systems.

For convenience, voltage equations of the PMSG are expressed in dq-frame [1]:

$$\begin{aligned} \upsilon\_{\eta} &= -r\_{s}\dot{i}\_{q} - L\_{q}\dot{i}\_{q} - \alpha \gamma L\_{d}\dot{i}\_{d} + \alpha \gamma \lambda\_{m} \\ \upsilon\_{d} &= -r\_{s}\dot{i}\_{d} - L\_{d}\dot{i}\_{d} + \alpha \gamma L\_{q}\dot{i}\_{q} \end{aligned} \tag{1}$$

where *vd* and *vq* are dq-frame voltages, *id* and *iq* are dq-frame currents, *rs* is winding resistance, *Ld* and *Lq* are dq-frame inductances, ω*<sup>r</sup>* is electrical rotor speed and λ*<sup>m</sup>* is the flux linkage established by the permanent magnet.

A PMSG can produce electromagnetic torque as follows:

$$T\_{\varepsilon} = \left(\frac{3}{2}\right) \left(\frac{P}{2}\right) \left[\lambda\_m i\_q + (L\_q - L\_d) i\_q i\_d\right]. \tag{2}$$

The electromagnetic torque is related to electrical rotor speed in the mechanical equation:

$$T\_{\varepsilon} = -J\left(\frac{2}{P}\right)\dot{\omega}\_{\tau} - B\left(\frac{2}{P}\right)\omega\_{\tau} + T\_I \tag{3}$$

where *J* is the inertia of a PMSG, *B* is the damping coefficient of a PMSG, *TI* is the input torque and *P* is magnetic pole number.

If *id* = 0 is satisfied in Equation (2), then the state equations of a PMSG are obtained:

$$\begin{array}{l}\dot{i}\_{q} = \frac{1}{L\_{s}}(-\upsilon\_{q} - r\_{s}\dot{i}\_{q} - \alpha \nu\_{L}\dot{i}\_{d} + \alpha \nu\_{r}\Lambda\_{m})\\\dot{i}\_{d} = \frac{1}{L\_{s}}(-\upsilon\_{d} - r\_{s}\dot{i}\_{d} + \alpha \nu\_{L}\dot{i}\_{q})\\\dot{\omega}\_{r} = \frac{1}{T}\Big[\left(\frac{P}{2}\right)(T\_{I} - T\_{c}) - B\omega\_{r}\Big] = \left(-\frac{3P^{2}}{8T}\lambda\_{m}\dot{i}\_{q} - \frac{P}{2}\alpha\nu\_{r} + \frac{P}{2}T\_{I}\right)\end{array} \tag{4}$$

#### **3. Design of the IT2 T-S Fuzzy Current Controller**

Two new state variables are defined to guarantee current tracking capability:

$$\begin{array}{lclclcl} s\_1 = \int [r\_1 - i\_q] dt & & \dot{s}\_1 = r\_1 - i\_q \\ s\_2 = \int [r\_2 - i\_d] dt & & \dot{s}\_2 = r\_2 - i\_d \end{array} \tag{5}$$

where *r*<sup>1</sup> is the target value of *iq* and *r*<sup>2</sup> is target value of *id*.

Extended state equations of the PMSG are obtained by combing Equations (4) and (5):

$$
\begin{bmatrix}
\dot{\alpha}\_{r} \\
\dot{i}\_{q} \\
\dot{i}\_{q} \\
\dot{i}\_{d} \\
\dot{s}\_{1}
\end{bmatrix} = \begin{bmatrix}
\frac{L\_{q}}{L\_{d}}\dot{i}\_{q} & 0 & -\frac{r\_{s}}{L\_{d}} & 0 & 0 \\
0 & -1 & 0 & 0 & 0 \\
0 & 0 & -1 & 0 & 0
\end{bmatrix} \begin{bmatrix}
\alpha\_{r} \\
\dot{i}\_{q} \\
\dot{i}\_{d} \\
\dot{s}\_{1} \\
\dot{s}\_{2}
\end{bmatrix} + \begin{bmatrix}
0 & 0 \\
0 & -\frac{1}{L\_{d}} & 0 \\
0 & 0 \\
0 & 0
\end{bmatrix} \begin{bmatrix}
\upsilon\_{d} \\
0 \\
\upsilon\_{d}
\end{bmatrix} + \begin{bmatrix}
\frac{PT\_{l}}{2f} \\
0 \\
0 \\
r\_{1} \\
r\_{2}
\end{bmatrix} \tag{6}
$$

$$\implies \dot{\mathbf{x}}(t) = A\mathbf{x}(t) + Bu(t) + E\boldsymbol{v}(t) \tag{7}$$

where *v*(*t*) represent disturbances in PMSG systems.

The required output function is:

$$\mathbf{y}(t) = \begin{bmatrix} 0 & 1 & 0 & 0 & 0 \\ & 0 & 0 & 1 & 0 & 0 \end{bmatrix} \begin{bmatrix} \alpha\_r \\ i\_q \\ i\_d \\ s\_1 \\ s\_2 \end{bmatrix} = \mathbf{C} \mathbf{x}(t). \tag{8}$$

The IT2 T-S fuzzy models and IT2 T-S fuzzy controller are combined to form the IT2 T-S fuzzy control systems. The characteristics of IT2 T-S fuzzy control systems employ the upper membership function and lower membership function to represent the model uncertainty of a nonlinear system. In the nonlinear system matrix of the PMSG, state variables *iq* and ω*<sup>r</sup>* are found. Therefore, *iq* and ω*<sup>r</sup>* are selected as Antecedent *z*<sup>1</sup> and *z*2, respectively. The membership function of the IT2 T-S fuzzy is shown

in Figure 2. *zp* is the Antecedent variable. *ap*, *bp*, *cp* and *dp* are the boundaries of the upper and lower membership functions. In the developed PMSG system, the rating of *iq* is 18A and the rating of ω*<sup>r</sup>* is 754 rad/sec. In order to represent the model uncertainty of the PMSG system, 5% variation of the Antecedent variable is selected. Therefore, the Antecedent variables and boundaries of the PMSG are summarized in Table 1.

**Figure 2.** Membership function of the IT2 T-S fuzzy control systems.

**Table 1.** Antecedent variables and boundaries of the PMSG.


In the design of the IT2 T-S fuzzy current controller, the nonlinear PMSG system is represented via linear sub-systems according to the model rules of the IT2 T-S fuzzy models:

$$\begin{array}{l}\text{Model rules } i:\\\text{If } z\_1(t) \text{ is } \widetilde{M}\_1^i \text{ and } \cdots \text{ and } z\_p(t) \text{ is } \widetilde{M}\_{p\prime}^i\\\text{then } \dot{\mathbf{x}}(t) = A\_i \mathbf{x}(t) + B\_i \boldsymbol{\mu}(t) + E\_i \mathbf{v}(t),\\\mathbf{y} = \mathbb{C}\_i \mathbf{x}(t), \ i = 1, 2, \cdots, r\end{array} \tag{9}$$

where *<sup>M</sup>*\**<sup>i</sup> <sup>p</sup>* are IT2 fuzzy sets, *x*(*t*) are state variables, *u*(*t*) are control inputs, *Ai*, *Bi* are state and input matrices of the sub-system, *Ei* is a constant matrix and *r* = 4 is the number of rules. The firing strength of the *i*-th rule is represented as follows:

$$
\overleftarrow{w}^{i}(z(t)) = \left[\underline{w}\_{i}(z(t)), \overline{w}\_{i}(z(t))\right], i = 1, 2, \cdots, r \tag{10}
$$

and

$$
\overline{w}\_i(z(t)) = \prod\_{j=1}^p \overline{\mu}\_{\widetilde{M}\_j^i}(z\_j(t)) \tag{11}
$$

$$\underline{w}\_{i}(z(t)) = \prod\_{j=1}^{p} \underline{\mu}\_{\widetilde{M}\_{j}^{i}}(z\_{j}(t)) \tag{12}$$

$$
\overline{w}\_i(z(t)) \ge \underline{w}\_i(z(t)) \ge 0, \forall i \tag{13}
$$

where *wi*(*z*(*t*)) is the upper grade of membership, *wi* (*z*(*t*)) ≥ 0 is the lower grade of membership, μ*M*\**i j* (*zj*(*t*)) is the upper membership function and μ *M*\**i j* (*zj*(*t*)) is the lower membership function. Then, the inferred IT2 T-S fuzzy model can be described as:

$$\mathbf{x}(t) = m \frac{\sum\_{i=1}^{r} \overline{w}\_{i}(z(t)) (A\_{i}\mathbf{x}(t) + B\_{i}u(t) + E\_{i}v(t))}{\sum\_{i=1}^{r} \overline{w}\_{i}(z(t))} + n \frac{\sum\_{i=1}^{r} \underline{w}\_{i}(z(t)) (A\_{i}\mathbf{x}(t) + B\_{i}u(t) + E\_{i}v(t))}{\sum\_{i=1}^{r} \underline{w}\_{i}(z(t))} \tag{14}$$

$$y(t) = m \frac{\sum\_{i=1}^{r} \overline{w}\_i(z(t)) \langle \mathcal{C}x(t) \rangle}{\sum\_{i=1}^{r} \overline{w}\_i(z(t))} + n \frac{\sum\_{i=1}^{r} \underline{w}\_i(z(t)) \langle \mathcal{C}x(t) \rangle}{\sum\_{i=1}^{r} \underline{w}\_i(z(t))} \tag{15}$$

where *m* and *n* are tuning parameters.

The parallel distributed compensation (PDC) of the IT2 T-S fuzzy controllers corresponding to the model rules are:

*Control rules i* :*Ifz*1(*t*) *is <sup>M</sup>*\**<sup>i</sup>* <sup>1</sup>*and*...*and zp*(*t*)*is <sup>M</sup>*\**<sup>i</sup> <sup>p</sup>*,*then u*(*t*) = *Kix*(*t*), *i* = 1, 2, ··· ,*r* (16)

where *Ki* is the controller gain. The inferred IT2 T-S fuzzy controller can be expressed as:

$$u(t) = m \frac{\sum\_{i=1}^{r} \overline{w}\_i(z(t)) (\mathcal{K}x(t))}{\sum\_{i=1}^{r} \overline{w}\_i(z(t))} + n \frac{\sum\_{i=1}^{r} \underline{w}\_i(z(t)) (\mathcal{K}x(t))}{\sum\_{i=1}^{r} \underline{w}\_i(z(t))}. \tag{17}$$

The close loop IT2 T-S fuzzy control system can be obtained by substituting Equation (17) into Equation (14).

#### **4. Stability Analysis**

Define the *H*<sup>∞</sup> performance index:

$$\sup\_{\|\psi(t)\|\_{2}\neq 0} \frac{\|y(t)\|\_{2}}{\|\psi(t)\|\_{2}} \le \gamma \text{, } 0 \le \gamma \le 1\tag{18}$$

where γ represents disturbance suppression ability of the IT2 T-S fuzzy control system.

**Lemma** [26]: Assume there exists a positive definite matrix *X* ∈ *n*×*<sup>n</sup>* and matrix *Mi* ∈ *m*×*n*, which makes the following LMI condition feasible:

$$
\begin{bmatrix}
\Phi\_{\text{iii}} & \* & \* \\
(m+n)\mathbb{C}\_i X & 0 & I
\end{bmatrix} \ge 0,
$$

$$
\begin{bmatrix}
\Phi\_{\bar{i}\bar{j}\bar{j}} + \Phi\_{\bar{j}\bar{u}} & \* & \* \\
( (m+n)(\mathcal{C}\_i + \mathcal{C}\_{\bar{j}}) )X & 0 & 2I
\end{bmatrix} \succeq 0, i < j \tag{20}
$$

$$\begin{bmatrix} \Phi\_{ijk} + \Phi\_{ikj} & \* & \* \\ -(2mE\_i + n(E\_j + E\_k))^T & 2\gamma^2 I & 0 \\ (2m\mathbb{C}\_i + n(\mathbb{C}\_j + \mathbb{C}\_k))X & 0 & 2I \end{bmatrix} \succeq 0, j < k \tag{21}$$

where

$$X = P^{-1},\ M\_i = K\_i X \tag{22}$$

$$\Phi\_{\text{ijk}} = -m(A\_i X + X A\_i^T) - n(A\_j X + X A\_j^T) + m(B\_i M\_k + M\_k^T B\_i^T) + n(B\_j M\_k + M\_k^T B\_j^T). \tag{23}$$

Then, the designed current controller in Equation (17) will guarantee the current tracking capability, which means the current tracking error will converge to zero. The detailed proof procedure can be referred in [26]. The controller gain can be found through using *Ki* = *MiX*<sup>−</sup>1.

Let γ = 0.86, *m* = 0.4 and *n* = 0.5; the following controller gains are obtained from the LMI conditions listed in Equations (19) to (21):

$$\begin{aligned} K\_1 &= \begin{bmatrix} 1.47332 & -0.28632 & -0.7274 & -893.16 & -0.18312 \\ 0.0222 & -0.7851 & -0.315 & 0.002 & -72.6 \end{bmatrix} \\ K\_2 &= \begin{bmatrix} 1.3759 & -0.28656 & -0.2424 & -796.401 & -0.0618 \\ 0.0208 & -0.2617 & -0.31524 & 0.001 & -72.6 \end{bmatrix} \\ K\_3 &= \begin{bmatrix} 1.5066 & -0.28668 & -0.7274 & -1414.17 & -0.18312 \\ 0.0021 & -0.7851 & -0.31536 & 0.002 & -72.6 \end{bmatrix} \\ K\_4 &= \begin{bmatrix} 1.3759 & -0.2391 & -0.2425 & -796.401 & -0.0618 \\ 0.0007 & -0.2617 & -0.3156 & 0.001 & -69.3 \end{bmatrix} \end{aligned} \tag{24}$$

#### **5. Results and Discussions**

The parameters of the PMSG are listed in Table 2. Table 3 shows the specifications of the PMSG drive. Figure 3 demonstrates the experimental equipment of the PMSG system based on the IT2 T-S fuzzy systems. The adopted PMSG is the model YBL17B-200L manufactured by YELI electric & machinery Co., LTD, Taiwan. The oscilloscope is KEYSIGHT DSO-X 3014T (Keysight Technologies, Santa Rosa, CA, USA). The current measurement system includes current probe Tektronix TCP303 and amplifier Tektronix TCPA300. Constant and variable current commend experimental results are given as follows.

**Table 2.** Parameters of the PMSG.


**Table 3.** Specifications of the PMSG Drive.


**Figure 3.** Experimental equipment for the PMSG system based on the IT2 T-S fuzzy systems.

#### *5.1. Constant Current Command*

Let the current command be *i* ∗ *<sup>q</sup>* = 15 A, *i* ∗ *<sup>d</sup>* = 0 A. In this case, the three-phase winding current will be balanced with peak value 15 A. Figure 4a–c show the winding current waveforms measured at different generator speeds. Table 4 summarizes the error and total harmonic distortion (THD) of the measured waveforms. It can be found that the current command tracking error is less than 2.5% in all conditions. The winding current can exactly track the current command. Moreover, the THD is less than 2%, and this means winding currents are nearly sinusoidal.

**Figure 4.** Winding current waveforms by letting *i* ∗ *<sup>q</sup>* = 15 A and *i* ∗ *<sup>d</sup>* = 0 A at different generator speeds: (**a**) 600 rpm, (**b**) 900 rpm and (**c**) 1200 rpm.


**Table 4.** Error and THD of the measured waveforms.

To improve the scientific merit of this paper, the simulation is performed by PSIM Ver. 10. The same conditions with Figure 4a–c are simulated and shown in Figure 5a–c, respectively. From the simulation results, it can be found that the three-phase winding currents are balanced with peak value 15 A. Moreover, the THD is less than 2%; this means winding currents are nearly sinusoidal.

**Figure 5.** Simulation winding current waveforms by letting *i* ∗ *<sup>q</sup>* = 15 A and *i* ∗ *<sup>d</sup>* = 0 A at different generator speeds: (**a**) 600 rpm, (**b**) 900 rpm and (**c**) 1200 rpm.

#### *5.2. Variable Current Command*

First, the current command *i* ∗ *<sup>q</sup>* is stepped from 5 A to 10 A (*i* ∗ *<sup>d</sup>* = 0) at 900 rpm. Figure 6a is the step occurring at 60◦ and Figure 6b is the step occurring at 120◦. The results using the IT1 T-S fuzzy system are compared in Figure 6c,d, respectively. The overshoot and settling time of the IT1 and IT2 T-S fuzzy systems are summarized in Table 5. It can be found that the winding currents can achieve the step current command faster by using the proposed IT2 T-S fuzzy control. The overshoot of the IT2 T-S fuzzy system is also less than the IT1 T-S fuzzy system, which makes the developed system more reliable and flexible.


**Table 5.** Overshoot and settling time of the IT1 and IT2 T-S fuzzy systems.

Note: overshoot = (maximum current − desired current) settling time: the required time when current tracking error is less than 5%.

IT2 T-S fuzzy 2.03 A 1.03 ms

Second, the speed is 1200 rpm and the current command is going from zero to a constant value. Figure 7a,b show the winding current waveforms of the current command that goes from zero to a constant value. The current tracking capability is verified.

**Figure 6.** Winding current waveforms of the current command stepped from 5 A to 10 A at 900 rpm: (**a**) step occurs at 60◦ with IT2 T-S fuzzy, (**b**) step occurs at 120◦ with IT2 T-S fuzzy, (**c**) step occurs at 60◦ with IT1 T-S fuzzy and (**d**) step occurs at 120◦ with IT1 T-S fuzzy.

**Figure 7.** Winding current waveforms of the current command goes from zero to a constant value at 1200 rpm: (**a**) 10 A; (**b**) 15 A.

Finally, the speed is 900 rpm and 1500 rpm, respectively, and the current command is varying as 0 A→5 A→10 A→15 A→10 A→5 A. Figure 8a,b show the winding current waveforms, respectively. It is obvious that the winding current can track the variable current command very well. The output power of the PMSG can be adjusted by changing the current command in this situation.

**Figure 8.** Winding current waveforms of the current command is varying as 0 A→5 A→10 A→15 A→10 A→5 A: (**a**) 900 rpm; (**b**) 1500 rpm.

#### *5.3. Calculation Time*

The calculation time of the IT2 T-S fuzzy and the IT1 T-S fuzzy systems is compared in Figure 9a,b. The calculation time of the IT2 T-S fuzzy is 32.6 μs and the calculation time of the IT1 T-S fuzzy is 30 μs. The IT2 T-S fuzzy system requires 2.6 μs more than the IT1 T-S fuzzy system to process the control algorithm. The switching period of the developed system is 50 μs. Therefore, the proposed algorithm is acceptable in implementation.

**Figure 9.** Calculation time of the IT2 T-S fuzzy and the IT1 T-S fuzzy systems: (**a**) IT2 T-S fuzzy; (**b**) IT1 T-S fuzzy.

#### **6. Conclusions**

The current control of the PMSG was designed and implemented based on the IT2 T-S fuzzy systems. First, the system configuration and the dynamic model were introduced. Next, the current

controller was designed based on the IT2 T-S fuzzy models. The IT2 T-S fuzzy control system was implemented to consider the uncertainty of nonlinear systems. The stability analysis and detailed design process were also demonstrated. The controller gain could be found by using the LMI conditions. Furthermore, the experimental equipment of the PMSG was illustrated. Experimental results, including constant current command, variable current command and calculation time, were demonstrated to verify the performance of the proposed current control. Simulation results were also performed in constant current command to improve the scientific merit of this paper. The IT2 T-S fuzzy system is more complex than the IT1 T-S fuzzy system. However, its calculation time is acceptable. Furthermore, its overshoot and settling time under the current command variation is better than the IT1 T-S fuzzy system. The output power of a PMSG can be adjusted by changing the peak value of three-phase balanced winding currents. In the future, the sensorless control combined with the IT2 T-S fuzzy systems may be adopted to increase the practicability of the PMSG system.

**Author Contributions:** The authors contributed equally to this work.

**Funding:** This research was supported by the Ministry of Science and Technology, Taiwan, under the Grant of MOST 108-2623-E-194-001-D.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **The Optimal Energy Dispatch of Cogeneration Systems in a Liberty Market**

**Whei-Min Lin 1, Chung-Yuen Yang 1, Chia-Sheng Tu 2, Hsi-Shan Huang <sup>2</sup> and Ming-Tang Tsai 3,\***


Received: 20 June 2019; Accepted: 20 July 2019; Published: 25 July 2019

**Abstract:** This paper proposes a novel approach toward solving the optimal energy dispatch of cogeneration systems under a liberty market in consideration of power transfer, cost of exhausted carbon, and the operation condition restrictions required to attain maximal profit. This paper investigates the cogeneration systems of industrial users and collects fuel consumption data and data concerning the steam output of boilers. On the basis of the relation between the fuel enthalpy and steam output, the Least Squares Support Vector Machine (LSSVM) is used to derive boiler and turbine Input/Output (I/O) operation models to provide fuel cost functions. The CO2 emission of pollutants generated by various types of units is also calculated. The objective function is formulated as a maximal profit model that includes profit from steam sold, profit from electricity sold, fuel costs, costs of exhausting carbon, wheeling costs, and water costs. By considering Time-of-Use (TOU) and carbon trading prices, the profits of a cogeneration system in different scenarios are evaluated. By integrating the Ant Colony Optimization (ACO) and Genetic Algorithm (GA), an Enhanced ACO (EACO) is proposed to come up with the most efficient model. The EACO uses a crossover and mutation mechanism to alleviate the local optimal solution problem, and to seek a system that offers an overall global solution using competition and selection procedures. Results show that these mechanisms provide a good direction for the energy trading operations of a cogeneration system. This approach also provides a better guide for operation dispatch to use in determining the benefits accounting for both cost and the environment in a liberty market.

**Keywords:** cogeneration systems; Time-of-Use (TOU); CO2 emission; Ant Colony Optimization

#### **1. Introduction**

Cogeneration systems, which combine heat and power (CHP) systems, have previously been extensively applied in industry. They offer an economic strategy providing both heat and power, which can then be passed on to buyers. Cogeneration systems offer a significant advantage when it comes to consideration of environmental issues. They are used as a distributed energy source, which can simultaneously sell both thermal steam and electricity to other industries. They can also be constructed in urban areas and used as distributed energy resources in microgrids [1–3]. In recent decades, consolidated cogeneration solutions have been used in industrial applications [4], while cogeneration system applications continue to grow. However, more experience is required in order to achieve the most efficient and energy-saving operation of these systems. To improve the competiveness of cogeneration systems in a liberalized market, an efficient tool for achieving the optimal operation of these systems must be developed.

To date, several efficiency strategies have been developed to achieve this optimal operation [5–19]. Ref. [5] presented a generalized network programming (GNP) to perform the economic dispatch of electricity and steam in a cogeneration plant. Ref. [6] presented an dispatching scheme which

economically transfers energy between facilities and utilities. An optimal operation of the cogeneration system is proposed, which will integrate energy into the electricity grid by using the decision-making technique [7]. Ref. [8] is used to assess the potential process of using micro-cogeneration systems based on Stirling engines. The results demonstrate that a numerical analysis of the Stirling engine can accurately indicate the operation of the actual machine. Ref. [9] used a non-linear programming with Time-of-Use (TOU) rates considered during operation. Ref. [10] presented the operation of steam turbines experiencing blades failures during peak load times of the summer months at a cogeneration plant. The author of [11] applied some possible technologies to integrate pulp and paper production within the context of a high-efficiency cogeneration system. Grey wolf optimization [12] and cuckoo search algorithms [13] are proposed to simultaneously solve the economic logistics of using a combined heat and power system. Ref. [14,15] presented the suggested economical operation of a cogeneration system under control, with resultant multi-pollutants from a fossil-fuels-based thermal energy generation. The author of [16] introduced an original framework based on identifying the characteristics of small-scale and large-scale uncertainties, whereby a comprehensive approach based on multiple time frames was formulated. Ref. [17] developed a tool for long-term optimization of cogeneration systems based on mixed integer linear-programming and a Lagrangian relation. Ref. [18] proposed an enhanced immune algorithm to solve the scheduling of cogeneration plants in a deregulated market. Ref. [19] addressed an optimal strategy for the daily energy exchange of a 22-MW combined-cycle cogeneration plant in a liberty market.

One of the key issues of a cogeneration operation is heat and power modeling. In the papers described above, pure power dispatch was a major objective. Inevitably, though, more design objectives coupled with higher constraints will have to be incorporated. The energy trading dispatch of cogeneration systems is a complicated process, especially when the solution is being sought in a world of uncertainty. Conventional methods have thus become more difficult to solve. Recently, artificial intelligence (AI) has been applied in the economic dispatch of cogeneration systems [20–23]. The strategies proposed by AI algorithms must consider computer execution efficiency and a large computing space. Conventional algorithms may be faster, but are very often limited by the problem structure, and may diverge, or lead to a local minimum. This paper therefore proposes an Enhanced Ant Colony Optimization (EACO) to solve the energy trading dispatch of cogeneration systems.

Ant Colony Optimization (ACO) applies the activity characteristics of biotic populations to search optimization problems [24,25]. When ants are foraging, they not only refer to their own information but also learn from the most efficient ants in order to correct their route. They learn and exchange their information to search for the shortest route between their colony and food sources, and pass this information on until the whole ant colony reaches optimal status. The advantages of the ACO algorithm are that individual solutions within a range of possible solutions can converge to discover the optimal solution through a small number of evolution iterations. ACO has previously been applied to the economic dispatch of power systems [26–30]; however, while ACO is good at global searches, the populations produced are still a dilemma. In this paper, an EACO algorithm is proposed to improve this search ability. In the EACO, the crossover and mutation mechanisms [31] are used to generate offspring equipped to escape the local optimum. This paper proposes the use of EACO to solve the energy trading dispatch of the cogeneration systems by considering the TOU rate [32]. The different carbon prices of CO2 emissions are also simulated and analyzed in the energy trading dispatch of these cogeneration systems. It can show the performance of the energy trading dispatch of the cogeneration systems to obtain the maximal profit.

#### **2. Problem Formulation**

Figure 1 shows M back pressure steam turbines, N extraction condensing steam turbines, and K high pressure steam boilers, where the high and medium pressure steam systems are connected by a common pipeline. Part of the generated electricity is supplied to the service power, and the excess

electricity is sold to the power company. Input/Output (I/O) models of the boiler and turbine are described as follows.

**Figure 1.** Energy flow of a cogeneration system.

#### *2.1. Operation Models of Boilers*

By using the Least-Square Support Vector Machine (LSSVM) [33], the operation models of boilers can be calculated from boiler operational records as shown in Figure 2. The temperature, pressure, fuel consumption, and steam generation for a real cogeneration system are measured from the operational data of boilers. The LSSVM is used for model training, and the input layer data are transferred to the output layer through the LSSVM. In general, using the Radial Basis Function Network (RBFN) kernel function, *K*(*x*, *y*) = *e*(−σ2|*x*−*y*<sup>|</sup> <sup>2</sup>), can yield a good prediction of the LSSVM model. Therefore, we adopt the RBFN kernel function as the kernel function of the LSSVM model. The error is calculated by using Mean Absolute Percentage Error (MAPE) as shown in Equation (1):

$$\text{MAPE} = \frac{1}{T} \sum\_{t=1}^{T} \frac{\left| S\_t^A - S\_t^F \right|}{S\_t^A} \times 100\% \tag{1}$$

where *SA <sup>t</sup>* is *<sup>t</sup>* actual operation data to be predicted *<sup>S</sup><sup>F</sup> <sup>t</sup>* is the *t* operation data constructed with LSSVM, and *T* is total training time.

**Figure 2.** Operation model of boilers using the LSSVM.

#### *2.2. Operation Model of Steam Turbines*

The steam from back-pressure steam turbines flows to the single inlet. After the high temperature and high pressure steam enters the steam turbine, the pressure reduces, the volume expands, and the temperature reduces. The steam in the outlet end is the 20 kg/300 ◦C medium pressure steam required for the process. The operation of the back-pressure steam turbine relies on the relationship between steam flow at the outlet and the generated electricity; the operation model of the back-pressure steam turbine is constructed by LSSVM, as shown in Figure 3.

**Figure 3.** Operation model of the back-pressure steam turbine using the LSSVM.

Extraction condensing turbines are different from the backpressure turbines, meaning the extraction condensing turbines have a single steam inlet. In the multi sections, after extraction of the middle/low pressure steam and exhaust of the steam at the final section, the condensing turbines are shown as in Figure 1. LSSVM is used to construct the generated electricity functions between the process steam outlet flow and the steam flow at the condensing section. The operation model constructed by LSSVM is shown in Figure 4.

**Figure 4.** Operation model of the extraction condensing steam turbines using LSSVM.

#### *2.3. Emission Model*

*CO*<sup>2</sup> emission models may be defined based upon the amount of fuel consumed. The model of emission for *CO*<sup>2</sup> is formulated by the IPCC [34] as:

$$E\_{co2,i}(t) = \mathcal{H}\Big(P\_{\mathcal{S}^i}(t)\Big) \times 4.1868 \times \frac{44}{12} \times \text{CEP}\_i \times \text{COR}\_i \tag{2}$$

H - *Pg*,*i*(*t*) = *di* + *eiPgi*(*t*) + *fiP*<sup>2</sup> *gi*(*t*) <sup>+</sup> *giP*<sup>3</sup> *gi*(*t*), approximated by the three order function gives the thermal conductivity of each kind of unit. *di*, *ei*, *fi*, *gi* are the coefficients of the emission of unit *i*. *CEPi* is the carbon emission parameter of unit *i* (21.1 kg-C/GJ for oil, 25.8 kg-C/GJ for coal, 15.3 kg-C/GJ for natural gas) and *CORi* is the carbon oxidizing rate of unit *i* (0.99 for oil, 0.98 for coal, 0.995 for natural gas).

#### *2.4. Objective Function and Constraints*

The purpose of the proposed scheme is to maximize profit while satisfying operational constraints. The objective function, including profit from steam sold, profit from electricity sold, fuel costs, emissions costs, wheeling costs, and water costs, is formulated in Equation (3):

$$\max T\mathcal{C}(\cdot) = \sum\_{t=1}^{T} \left[ T\mathcal{C}\_{\text{stam}}(t) + T\mathcal{C}\_{\text{clart}}(t) - T\mathcal{C}\_{\text{fucl}}(t) - T\mathcal{C}\_{\text{emiss}}(t) - T\mathcal{C}\_{\text{tram}}(t) - T\mathcal{C}\_{\text{nutz}}(t) \right] \tag{3}$$

(1) The profit from thermal steam sold:

$$TC\_{steam}(t) = S\_{out}(t) \times Stam\_{cost} \tag{4}$$

(2) The profit from electricity sold:

$$TC\_{\text{elect}}(t) = \text{TOU}(t) \times P\_{\text{tie}}(t) \tag{5}$$

(3) The cost of fuel:

$$TC\_{fuel}(t) = \sum\_{i=1}^{K} F\_{bi}(\mathcal{S}\_{bi}(t)) \times \text{Fuel\\_cos}\, t \tag{6}$$

(4) Emission costs:

$$TC\_{\text{emiss}}(t) = \mathbb{C}\_{\mathbb{C}} \times \sum\_{i=1}^{K} E\_{\text{CO2},i}(t) \tag{7}$$

(5) Wheeling costs:

$$TC\_{tmn}(t) = \mathbf{WJC} \times P\_{tic}(t) \tag{8}$$

(6) The cost of pure water:

$$TC\_{\text{water}}(t) = \mathsf{WC} \times \mathsf{W}\_{b}(t) \tag{9}$$

Cc: The charged emission price for CO2.(NT\$235/ton) [35]; *Fbi*(*Sbi*(*t*)): consumed enthalpy of the *i* steam boiler at *t* hour; *Fuel\_cost*: The fuel cost(NT\$/MBTU) for coal, gas, and oil; *Sout*(*t*): The thermal steam sold at time *t* (ton/h); *Steamcost*: The price of thermal steam sold (NT\$/ton); *Ptie*(*t*): The electricity sold to utility at time *t*; *TOU*(*t*): Time-of-Use rate (NT\$/KWH), as shown in Table 1 [32]; *Wb*(*t*): The water used by boilers at time *t* (ton/h); *WC*: The cost of water (NT\$/ton); WUC: the wheeling cost (NT\$/MWH).


**Table 1.** The time-of-use rates of Taipower Company.

Level 1: power exported at under 20% rated capacity; Level 2: power exported at over 20% rated capacity.

The operation constraints are considered as follows:

(1) Power balance in the power system:

$$\sum\_{i=1}^{M+N} P\_{\mathcal{g}^i}(t) - P\_{it\epsilon}(t) - P\_{load}(t) = 0 \tag{10}$$

(2) Steam balance for boilers, turbines, sold, and industrial processes:

$$\sum\_{i=1}^{K} S\_{hi}(t) - D\_h(t) - \sum\_{i=1}^{M+N} S\_{lii}(t) = 0 \tag{11}$$

$$\sum\_{i=1}^{M+N} S\_{hi}(t) - \sum\_{i=1}^{M+N} S\_{mi}(t) - \sum\_{i=1}^{N} S\_{wi}(t) = 0 \tag{12}$$

$$\sum\_{i=1}^{M+N} S\_{mi}(t) - D\_m(t) - S\_{out}(t) = 0\tag{13}$$

(3) Operation constraints for boilers, steam turbines and power generation:

$$S\_{bi}^{\min} \le S\_{bi}(t) \le S\_{bi}^{\max} \quad , \ i = 1, 2, 3, \dots, K \tag{14}$$

$$S\_{hi}^{\min} \le S\_{hi}(t) \le S\_{hi}^{\max} \qquad , \text{ } i = 1, 2, 3, \dots, M+N \tag{15}$$

$$S\_{\rm mi}^{\rm min} \le S\_{\rm mi}(t) \le S\_{\rm mi}^{\rm max} \quad , \ t = 1, 2, 3, \dots, M+N \tag{16}$$

$$S\_{wi}^{min} \le S\_{wi}(t) \le S\_{wi}^{max} \quad , \ i = 1, 2, 3, \dots, N \tag{17}$$

$$P\_{\mathcal{S}^i}^{\min} \le P\_{\mathcal{S}^i}(t) \le P\_{\mathcal{S}^i}^{\max} \quad , \ i = 1, 2, 3, \dots, M+N \tag{18}$$

*Dh*(*t*), *Dm*(*t*): The high/medium pressure steam demands of industry at time t (T/h); *Pload*(*t*): The load of cogeneration system at time t; *S*min *bi* , *<sup>S</sup>*max *bi* : Minimal/maximal limits of steam for boiler *<sup>i</sup>*; *<sup>S</sup>*min *hi* , *S*max *hi* : Minimal/maximal limits of high pressure steam for turbine *<sup>i</sup>*; *<sup>S</sup>*min *mi* , *<sup>S</sup>*max *mi* : Minimal/maximal limits of medium pressure extraction steam for turbine *i*; *S*min *wi* , *<sup>S</sup>*max *wi* : Minimal/maximal limits of medium pressure exhausted steam for turbine *i*; *P*min *gi* , *<sup>P</sup>*max *gi* : Minimal/maximal limits of the generated power for turbine *i*.

#### **3. Solution Algorithm**

This study proposes an EACO, which combines the ACO and Genetic Algorithm (GA), in order to achieve the optimal energy trading dispatch of a cogeneration system. Crossover and mutation mechanisms are integrated into the ACO procedure, and serve to generate offspring in order to escape from the local optimum. The EACO procedure applied in the energy trading dispatch of a cogeneration system is described as follows.

#### **(1) Input Data**

Input data includes high/medium steam demand, internal load, plant type, plant capacity, and number of plants.

#### **(2) Set EACO Parameters**

EACO parameters include the population of ants (k), the number of generations (G), initial pheromone (τ<sup>0</sup> = 0.1), the relative influence of the pheromone trail (α = 1), the relative influence of the heuristic information (β = 2), and the pheromone evaporation rate (ρ = 0.5).

#### **(3) Initialized Individuals**

R*i <sup>s</sup>* <sup>=</sup> , *Si b* , *Si <sup>m</sup>*, *S<sup>i</sup> w* is an individual, *i* = 1, 2, ... , *k*. *k*, which is the number of ants, is set to 30 in our study. *s* is the number of parameters. All individuals are set between the minimal and maximal limits with a uniform distribution as shown in Equation (19). The fitness score of each R*<sup>i</sup> <sup>s</sup>* is obtained by calculating the objective function (TC(·)) by considering Equations (10)~(18). The fitness values were arranged in descending order from the maximum (*TC*(R*<sup>i</sup> s*)*max*) to the minimum (*TC*(R*<sup>i</sup> <sup>s</sup>*)*min*).

$$\mathcal{R}\_{\rm s}^{i} = \mathcal{R}\_{\rm s\,min}^{i} + \text{Rand} \times \left(\mathcal{R}\_{\rm s\,max}^{i} - \mathcal{R}\_{\rm s\,min}^{i}\right) \tag{19}$$

*Rand*: The uniform random number in (0,1).

#### **(4) The State Transition Rule**

The state-based ants are generated according to the level of pheromone and constrained conditions. The transition probability for the *k* − *ant* from one state *s* to the next *j* is at the *t* − *th* interval, as given in Equation (20):

$$P\_{t,sf}^{k}(\mathbf{g}) = \begin{cases} \frac{\left[\tau\_{t,j}(\mathbf{g})\right]^a \times \left[\eta\_{t,j}(\mathbf{g})\right]^\beta}{\Sigma\_{lk\mathbf{N}\_l^k(s)} \left[\tau\_{t,k}(\mathbf{g})\right]^a \times \left[\eta\_{t,l}(\mathbf{g})\right]^\beta}, & \text{if } j \in \mathbf{N}\_l^k(s) \\\\ 0 & \text{, otherwise} \end{cases} \tag{20}$$

where η*t*,*sj*(*g*) and η*t*,*sl*(*g*) are the inverse of the edge distance at the *g* − *th* generation, which are defined as Equations (21) and (22):

$$\eta\_{t;sl}(\mathcal{g}) = \frac{1}{\left| TC(\mathcal{R}\_{t,s}) - TC(\mathcal{R}\_{t,\text{optimal}}) \right|} \tag{21}$$

$$\eta\_{l;sl}(\mathcal{g}) = \frac{1}{\left| \text{TC}(\mathcal{R}\_{t,s}) - \text{TC}(\mathcal{R}\_{t,l}) \right|} \quad , \ l \in N\_l^k(s) \tag{22}$$

*TC*(·) is the objective function as given in Equation (3). *TC*- R*t*,*s* and *TC*- R*t*,*<sup>l</sup>* are the score of the *<sup>s</sup>* <sup>−</sup> *th* and *<sup>l</sup>* <sup>−</sup> *th* individuals at the *<sup>t</sup>* <sup>−</sup> *th* interval, and *TC* <sup>R</sup>*t*,*optimal* is the optimal fitness score at the *t* − *th* interval. *N<sup>k</sup> <sup>t</sup>*(*s*) is the number of feasible individuals at the *t* − *th* interval.

τ*t*,*sj*(*g*) and τ*t*,*sl*(*g*) are the pheromone intensity on edge *(s, j)* and edge *(s, l)* at the *g* − *th* generation. Ant k positioned on state *s* chooses to move to the next state by taking account of τ*t*,*sl* and η*t,sl*. When the value of τ*t,sl* increases, this indicates there has been a lot of traffic on this path and it is therefore more desirable in order to reach the optimal solution. When the value of η*t*,*sl* increases, it indicates that the current state should have a higher probability. Each stage contains several states, while the order of states selected at each stage can be combined as an achievable path deemed to be a feasible solution to the problem.

#### **(5) Ant Reproduction**

New ants are generated by the scheme of crossover and mutation. Crossover is a structured recombination operation that exchanges two individual ants. Mutation is the occasional random alteration of an individual. The crossover and mutation scheme is described as follows:

	- (a) If *rang* > *CV*(*g*) : Mutation is used;
	- (b) If *rang* > *CV*(*g*) : Crossover is used. *Rand*: the uniform random in [0, 1]; *CV*(*g*): the control variable between 0.1 to 0.9. The initial value set to 0.5; *g*: the current iteration number.

$$CV(\text{g} + 1) = CV(\text{g}) - RP \tag{23}$$

where *RP* <sup>=</sup> <sup>|</sup>*CV*(*g*)−*CV*(*g*−1)<sup>|</sup> <sup>2</sup> is the regulated parameter. When *RP* is added in the crossover process, the higher probability for crossover operation will produce the next offspring.

(iii) If *TC*- *t*,*optimal* comes from mutation, the control parameter *CV*(*g* + 1) will be increased as shown in Equation (24):

$$CV(\text{g} + 1) = CV(\text{g}) + D \tag{24}$$

Similarly, when *RP* is added in the mutation process, the higher probability for mutation operation will produce the next offspring.

(iv) If *TCmin*(*g* − 1) = *TCmin*(*g*), the control variable needs to hold back. If *CV*(*g*) > *CV*(*g* − 1), we have:

$$CV(\text{g} + 1) = CV(\text{g}) - D \tag{25}$$

otherwise,

$$CV(\emptyset + 1) = CV(\emptyset) + D \tag{26}$$

The crossover operator proceeds to exchange two individual ants by random. R*t*,*<sup>s</sup>* and R*t*,*l*, which are the *s* − *th* and *l* − *th* individual ants at the *t* − *th* interval, are exchanged by the crossover operator. The mutation operation is carried out to produce another individual ant. Each individual ant is mutated and created to a new individual ant by using (27).

$$\mathcal{R}\_{t,s}^{j+1} = \mathcal{R}\_{t,s}^{j} + \mathcal{N}(0, \sigma^2) \tag{27}$$

*N* - 0, σ<sup>2</sup> represents a Gaussian random variable with mean 0 and variable σ2. σ<sup>2</sup> can be calculated by:

$$
\sigma = \beta \times \left( \mathcal{R}\_{t, \text{s. max}}^{\text{j}} - \mathcal{R}\_{t, \text{s. min}}^{\text{j}} \right) \times \frac{T \mathbb{C} \left( \cdot \right)^{\text{k}}}{T \mathbb{C} \left( \cdot \right)\_{\text{max}}^{\text{k}}} \tag{28}
$$

β, which is a mutation factor at the *j*-th generation, is set within [0, 1].

#### **(6) Update the Pheromone**

While building a solution to the problem, the pheromone of a visited route can be dynamically adjusted by Equation (29). This process is called the "local pheromone-updating rule":

$$
\pi\_{t,sj}^{k+1} = (1 - \rho)\tau\_{t,sj}^k + \Delta \tau\_{t,sj}^k \tag{29}
$$

ρ is the constant of pheromone intensity (0 ≤ ρ ≤ 1) and Δτ*<sup>k</sup> <sup>t</sup>*,*sj* is the deviation of pheromone intensity on edge (*s, j*) at the *t* − *th* interval, as shown in Equation (30):

$$
\Delta \tau\_{t, \text{st}}^k = \left\{ \begin{array}{ll} \mathbb{Q}/c\_{t, \text{sj}}^k & , \quad \text{the } \text{path}(\text{s}, j) \text{ for } k-\text{th} \text{ and} \\\\ 0 & , \quad \text{other} \end{array} \right\} \tag{30}
$$

*Q* is the release rate of pheromone (0 ≤ *Q* ≤ 1) and Δ*e<sup>k</sup> <sup>t</sup>*,*sj* is the path error (*s, j*) for the *k-th* ant.

#### **(7) Stopping Rule**

If a pre-specified stopping condition is satisfied, the run must be stopped and the results outputted; otherwise, return to step 4. In this study, the stopping rule is set at 300 generations.

#### **4. Case Study**

The proposed algorithm was tested with three back-pressure steam turbines, four extraction condenser steam turbines and seven steam boilers using gas, oil, and coal as fuel. Steam generation was measured in the field. All facilities including generators, boilers, and steam turbines have their capacity limitation. The limits of all facilities are introduced in Tables 2 and 3.


**Table 2.** Maximal and minimal limits of boiler flows and generators.


**Table 3.** Steam output limits.

#### *4.1. Modeling Tests for Boilers and Steam Turbines*

In this paper, the operation data of the boilers in the cogeneration system are recorded during the periods of each normal working day. The number of operation data samples for each boiler is 60; these are used to establish the operational models of boilers 1–7. The error results are shown in Table 4, and the example for operation models of boilers 1–2 are shown in Figure 5.

**Table 4.** The error results for the operation model of boilers.


**Figure 5.** Example of the operation models of boilers 1–2.

In this paper, the operation data of the backpressure steam turbines (ST1~3) and extraction condensing steam turbines (ST4~7) in the cogeneration system are recorded and used by LS and LSSVM to establish the operation models of the steam turbines (ST1~7). The error results are shown in Table 5, and the example for the operation curves of steam turbines ST1 and ST6 are shown in Figure 6. The error results for the operation model of steam turbines is from 1.3916% to 2.1475%. It can be proved that the accuracy is reliable.

**Table 5.** The error results for the operation model of steam turbines.


**Figure 6.** Example for the operation models of steam turbines (ST1/ST6).

#### *4.2. Results at Di*ff*erent TOU Intervals*

EACO was used to solve the energy dispatch of the cogeneration systems under a liberalized market. In the study system, the internal load was 30 MW and the steam demands of the industry for high- and medium-pressure were 60 T/H and 600 T/H, respectively. Table 6 shows the simulation results at the different TOU intervals. Table 6 indicates that all constraints were met. Since the electricity price is higher at the peak period, the operation strategy of cogeneration systems sold more power to the utility in order to achieve greater benefit. During off-peak periods, the cogeneration systems sell more thermal steam power to attain further benefits. It can be shown that the operation dispatch of the cogeneration systems during the peak period generated 273.321 MW and sold about 243.321 MW to

the utility. During the peak period, all generators produced their highest outputs, while tending to sell the thermal steam to industries during semi-peak and off-peak periods. The TOU rate plays an important role in the economic dispatch of a cogeneration system.


**Table 6.** Dispatch results at different TOU intervals.

Table 7 gives a profit analysis of cogeneration systems at different TOU intervals. The profits during peak period, semi-peak period, and off-peak peroid are 267,678.05, 13,606.32, and −121,047.33, respectively. The profit is lost during the off-peak period. As the cogeneration systems sell more power to utilities during the peak period, greater profits are realized. In the off-peak period, cogeneration systems sell more thermal steam to industries, thereby optimizing the energy trading dispatch. From Table 7, it is noted that the TOU rate influences the profits of the cogeneration system.

**Table 7.** Profit analysis of cogeneration systems at different TOU intervals. Unit: NT\$/H.


#### *4.3. Convergence Test*

Table 8 shows the comparisons of EP, GA, PSO, ACO, and EACO during different periods. An IBM PC with a P-IV2.0 GHz CPU and 512 MB SDRAM was used for this test. From this, the improvement of the EACO over other algorithms is clear. Figure 7 shows the convergent characteristics of EP, GA, PSO, ACO, and EACO during the peak period. The average execution times for EACO and ACO were only 1.85 s and 1.67 s, respectively. Although the executed performance of EACO was subtle, it showed the capacity of EACO to explore a solution more likely to offer maximum benefit.

**Table 8.** Comparison of EP, GA, PSO, ACO, and EACO algorithms. Unit: NT\$/H.


**Figure 7.** Convergence characteristics of EP, GA, PSO, ACO, and EACO during the peak period.

#### *4.4. Robustness Test*

Four algorithms including EP, GA, PSO, and ACO, were also tested at the peak period; the results are shown in Table 9. Each algorithm was executed using 100 trials during the robustness test. The number of optimum and average converged profit for EACO are 87 and NT\$261,618.29. It is seen that EACO offers a greater ability to attain maximum profits and a higher probability of finding the best solution.


**Table 9.** Robustness test for EP, GA, PSO, ACO, and EACO algorithms in a peak period.

#### *4.5. The Influence of Carbon Price*

Table 10 shows the impact of various carbon prices experienced during the peak period. Table 10 suggests that if carbon prices are higher, CO2 emissions will decrease. Similarly, due to the fuel type for extraction-condenser turbines being oil or coal, power generation will be reduced depending upon the carbon price. The purpose of various carbon prices here is to illustrate the tradoff between profit and emission costs, and also to show that generators may more economically dispatch trade electricity or CO2 emission to find better profit.

**Table 10.** The impact of various carbon prices during the peak period.


Figure 8 shows the profit-emission tradeoff curves during the peak period. Figure 8 provides diversified alternatives for decision makers, showing the effects of various carbon prices. Replaced with the maximal allowable emission as a constraint, an appropriate decision can be chosen to satisfy the desired level of profit and emission costs.

**Figure 8.** Profit-emission tradeoff curves during the peak period.

#### **5. Conclusions**

In this paper, an EACO is proposed to maximize the profit of energy trading using a cogeneration system. The objective function is formulated based upon a maximal profit model, which includes profit from steam sold, profit from electricity sold, fuel costs, CO2 emission costs, wheeling costs, and water costs. By considering the various carbon prices, the profits of the energy trading dispatches are evaluated while considering three different TOU scenarios. The effectiveness of the EACO is demonstrated and simulated on a real cogeneration system. Our analysis points to expectations of the TOU rate or carbon price for the energy trading dispatch. With the advantages of both heuristic ideals and ACO, EACO has threefold conventional ideals: the complicated problem is solvable, with a better performance than ACO, and the more likelihood to get a global optimum than heuristic methods. The results indicate that both provide good tools for determining the optimum energy trading operation of a cogeneration system. This shows that the tradeoff between investment cost and environmental protection can be clearly predetermined in the liberty market. EACO also has great potential to be further applied to many ill-conditioned problems in power system planning and operations.

**Author Contributions:** W.-M.L. generalized the novel algorithms and designed the system planning projects. C.-Y.Y. used system parameters and models for the simulation test. M.-T.T. performed the editing and experimental model results. H.-S.H. written the orginal draft perparation. C.-S.T. provided hardware tools and system model related materials. All the authors were involved in exploring system validation and results, and permitting the benefits of the published document.

**Funding:** This research was funded by Ministry of Science and Technology of R.O.C. grant number MOST 107-2221-E-230-008-MY3.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
