Multiobjective Genetic Algorithm-Based Optimization of PID Controller Parameters for Fuel Cell Voltage and Fuel Utilization

Abstract

Nowadays, given the great deal of fossil fuel consumption and associated environmental pollution, solid oxide fuel cells (SOFCs) have shown their great merits in terms of high energy conversion efficiency and low emissions as a stationary power source. To ensure power quality and efficiency, both the output voltage and fuel utilization of an SOFC should be tightly controlled. However, these two control objectives usually conflict with each other, making the controller design of an SOFC quite challenging and sophisticated. To this end, a multi-objective genetic algorithm (MOGA) was employed to tune the proportional–integral–derivative (PID) controller parameters through the following steps: (1) Identifying the SOFC system through a least squares method; (2) designing the control based on a relative gain array (RGA) analysis; and (3) applying the MOGA to a simulation to search for a set of optimal solutions. By comparing the control performance of the Pareto solutions, satisfactory control parameters were determined. The simulation results demonstrated that the proposed method could reduce the impact of disturbances and regulate output voltage and fuel utilization simultaneously (with strong robustness).

1. Introduction

Nowadays, with a great deal of fossil fuel consumption, many environmental issues have been raised. With the merits of low pollution and high efficiency [1], hydrogen energy-based fuel cells are considered to be ideal alternative energy sources, and they are widely used to supply power to stand-alone power systems and vehicles [2,3,4]. Among all fuel cells, solid oxide fuel cells (SOFCs) attract much attention for their multifuel capabilities and potential for cogeneration [5], as well as their ease of coordination with refinery hydrogen networks [6]. Numerous studies about SOFCs have been done. Rao et al. [7] analyzed the thermodynamic characteristics and techno-economic performance of a system consisting of an SOFC–gas turbine system, Xia et al. [8] studied low-temperature operations of SOFCs, and Wang et al. [9] modified the anode material to solve the problem of coke formation. Besides these, studies on efficient control strategies to ensure the safety, stability, and efficiency of SOFCs are also very important. In order to maintain the efficiency of an SOFC, fuel utilization in an SOFC must be kept between 70% and 90% [10]. In addition, to ensure the power quality of an SOFC, the output voltage must be controlled tightly. However, due to the nonlinearity and couplings of multiple variables and frequent load disturbances, the control of an SOFC is extremely hard, and the controller must be appropriately designed. As a result, designing an efficient control scheme is of great significance in the research on SOFCs.

To solve the problem of maintaining the output voltage and fuel utilization of SOFCs, many advanced control strategies have been proposed, such as model predictive control (MPC) [11,12,13], fuzzy proportional–integral–derivative (PID) control [14], fuzzy logic control [15], and neural network control [16]. All of these control strategies have been shown to obtain excellent control performance in numerous simulation studies: However, due to their computational complexity, these control strategies are difficult to implement in practical applications.

With the merits of a simple structure and high reliability, PID controllers are widely used in the process industry: However, to ensure control performance, the control parameters of PID controllers must be tuned well. To this end, numerous artificial intelligence (AI)-based optimization methods have been applied to tune PID controllers. For example, Qin, Sun, and Hua [17] used a fruit fly optimization algorithm (FOA) to optimize the PID controllers of a refrigeration system, Zamani et al. [18] applied a particle swarm optimization (PSO) algorithm to optimize the PID controller of an automatic voltage regulator (AVR), and Elbayomy et al. [19] optimized the PID controller of an electrohydraulic servo actuator system (EHSAS) to control the movable surface of a space vehicle. However, these references all used a single-objective optimization method, where the authors judged a controller’s performance through one objective function that weighted the sum of different performance indices, such as the integral absolute error (IAE), overshoot, and settling time. Although in the end this method achieved good results, it ignored the problem that the control objectives were usually conflicting and that the weight factors were decided subjectively by the author. Instead of using only one objective function, a Pareto-based multiobjective optimization algorithm involves two or more objective functions and is thus able to deal with the conflicting objectives of a control system without designating weight factors: Many researchers choose this method to solve optimization problems [20,21,22,23].

This paper aims to (1) identify a dynamic SOFC model as a linearized model, (2) apply two PID controllers to maintain the output voltage and fuel utilization of an SOFC under load disturbance, and (3) enhance the control system performance by optimizing the controller with multiobjective optimization. The remaining part of this paper is organized as follows: In Section 2, a dynamic SOFC model is introduced, while in Section 3, a linearized model is identified and a controller for the model is designed. Section 4 optimizes the controller parameters with genetic algorithm (GA)-based multiobjective optimization, and Section 5 applies the optimal controllers to the nonlinear system and analyzes the optimization results. Conclusions are drawn in Section 6.

2. Problem Description

2.1. Dynamic Model for an SOFC

The model presented in this paper was based on three submodels: An electrochemical model, a mass balance model, and an energy balance model. To simplify the model, six assumptions were made as follows:

• Both hydrogen and air are preheated to a specific temperature before entering the anode and cathode;
• Both hydrogen and air are in full contact with the anode and cathode, and the stoichiometric quantity of oxygen at the cathode is sufficiently large;
• A large stoichiometric quantity of oxygen exists at the cathode;
• The oxygen concentration in the air is 21%;
• The whole stack can be presented as a combination of individual stacks;
• Ideal gas laws are employed for both the fuel flow and air flow.

A structural diagram of this model is shown in Figure 1, where the fuel flow, air flow, oxygen ion flow, and electron flow are represented by red, blue, yellow, and green arrows, respectively.

2.1.1. Electrochemical Model

As we can tell from Figure 1, the SOFC consists of an anode and a cathode, and they are connected by a solid ceramic electrolyte. Through the electrolyte, oxygen ions are transported from the cathode to the anode, and electrons are transported from the anode to the cathode through the external electric circuit. The reaction on the anode and the cathode and the total reaction are shown as follows:

$$2H_2+2O^{2-}\to 2H_{2}O+4e^{-},$$

(1)

$$O_2+4e^{-}\to 2O^{2-},$$

(2)

$$2H_2+O_2\to H_{2}O.$$

(3)

Based on Nernst’s equation [14], the open circuit voltage of the fuel cell stack E_cell can be formulized as

$$E_{cell}=E_{0,cell}+\frac{R{T}_{cell}}{2F}\ln \left[\frac{p_{H_2}{p}_{O_2}^{0.5}}{p_{H_{2}O}}\right].$$

(4)

Here, E_0,cell can be calculated as

$$E_{0,cell}=E_{0,cell}^0-k_E\left(T-298\right),$$

(5)

where $E_{0,cell}^0$ is the reference voltage at 298 K and 1 atm of pressure.

The partial pressures in Equation (4), p_H2, p_O2, and p_H2O, are the average of the inlet and outlet partial pressures of the respective components, calculated as follows:

$$p_{H_2}=\frac{p_{H_2}^{in}+p_{H_2}^{out}}{2},$$

(6)

$$p_{O_2}=\frac{p_{O_2}^{in}+p_{O_2}^{out}}{2},$$

(7)

$$p_{H_{2}O}=\frac{p_{H_{2}O}^{in}+p_{H_{2}O}^{out}}{2}.$$

(8)

Due to ohmic loss, activation loss, and concentration loss, the output voltage of the fuel cell stack V_cell is lower than E_cell, which is calculated as

$$V_{cell}=E_{cell}-V_{ohm}-V_{act}-V_{con}.$$

(9)

The expression of the three voltage drops are given below:

$$V_{ohm}=I{A}_{cell}{R}_{ohm},$$

(10)

$$V_{act}=\frac{RT}{F}\left({\sinh }^{-1}\left(\frac{I}{2i_{0,an}}\right)+{\sinh }^{-1}\left(\frac{I}{2i_{0,ca}}\right)\right),$$

(11)

$$V_{con}=\frac{RT}{2F}\left(\ln \left(\frac{p_{H_{2}O,r}{p}_{H_2}}{p_{H_{2}O}{p}_{H_2,r}}\right)+\frac{1}{2}\ln \left(\frac{p_{O_2}}{p_{O_2,r}}\right)\right),$$

(12)

where R_ohm is the internal resistance of the fuel cell, and it is calculated as follows [24]:

$$R_{ohm}=\left(\frac{a_{ele}\exp \left(\frac{b_{ele}}{T}\right)}{A_{cell}}\right){\delta }_{ele}+\left(\frac{a_{itc}\exp \left(\frac{b_{itc}}{T}\right)}{A_{cell}}\right){\delta }_{itc}.$$

(13)

According to assumption 5, the output voltage of the fuel cell can be calculated with V_cell by

$$V_{out}=N{V}_{cell},$$

(14)

where N is the number of fuel cells in a series.

2.1.2. Mass Balance Model

According to assumption 5, the dynamic effective partial pressure of H₂, O₂, and H₂O can be presented as follows [25]:

$$\frac{V_{an}}{RT}\frac{d{p}_{H_2}^{ch}}{dt}=M_{H_2}^{in}-M_{H_2}^{out}-\frac{I}{2F},$$

(15)

$$\frac{V_{an}}{RT}\frac{d{p}_{H_{2}O}^{ch}}{dt}=M_{H_{2}O}^{in}-M_{H_{2}O}^{out}+\frac{I}{2F},$$

(16)

$$\frac{V_{ca}}{RT}\frac{d{p}_{O_2}^{ch}}{dt}=M_{O_2}^{in}-M_{O_2}^{out}-\frac{I}{2F},$$

(17)

while based on assumption 2, the inlet and outlet mass flow rate of H₂, O₂, and H₂O can be described as

$$\{\begin{array}{l}{M}_{H_2}^{in}=M_{an}{x}_{H_2}^{in}=M_{an}\frac{p_{H_2}^{in}}{p_{an}^{ch}}\\ M_{H_2}^{out}=M_{an}{x}_{H_2}^{out}=M_{an}\frac{p_{H_2}^{out}}{p_{an}^{ch}}\end{array},$$

(18)

$$\{\begin{array}{l}{M}_{H_{2}O}^{in}=M_{an}{x}_{H_{2}O}^{in}=M_{an}\frac{p_{H_{2}O}^{in}}{p_{an}^{ch}}\\ M_{H_{2}O}^{out}=M_{an}{x}_{H_{2}O}^{out}=M_{an}\frac{p_{H_{2}O}^{out}}{p_{an}^{ch}}\end{array},$$

(19)

$$\{\begin{array}{l}{M}_{O_2}^{in}=M_{ca}{x}_{O_2}^{in}=M_{an}\frac{p_{O_2}^{in}}{p_{ca}^{ch}}\\ M_{O_2}^{out}=M_{ca}{x}_{O_2}^{out}=M_{an}\frac{p_{O_2}^{out}}{p_{ca}^{ch}}\end{array}.$$

(20)

2.1.3. Energy Balance Model

As Figure 1 shows, hydrogen flows through the anode surface and changes into water, while air is provided by a central air supply tube (AST) and is forced to flow through the cathode. Due to the simplification of the model, the means of heat transfer inside the cell are radiation, convection, and mass flow [26].

Therefore, the energy balance equations for different parts of the cell are shown as follows:

• Cell tube:

  $$q_{cell,in}=q_{gen}=q_{chem}-q_{ele},$$

  (21)

  $$q_{cell,out}=q_{rad}+q_{conv,ann}+q_{flow,air,ann}+q_{conv,fuel}+q_{flow,fuel},$$

  (22)

  $$q_{cell,net}=q_{cell,in}-q_{cell,out}=m_{cell}{C}_{cell}\frac{d{T}_{cell}}{dt},$$

  (23)

  $$q_{rad}={\epsilon }_{AST}\sigma A_{AST,outer}\left(T_{cell}^4-T_{AST}^4\right),$$

  (24)

  $$q_{conv,ann}=h_{cell}{A}_{cell,inner}\left(T_{cell}-T_{air,ann}\right),$$

  (25)

  $$q_{flow,air,ann}=M_{air}{x}_{air}{C}_{air}\left(T_{air,ann,in}-T_{air,ann,out}\right),$$

  (26)

  $$q_{conv,fuel}=h_{cell}{A}_{cell,outer}\left(T_{cell}-T_{fuel}\right),$$

  (27)

  $$\begin{array}{ll}{q}_{flow,fuel}& =\left(M_{H_2,in}+M_{H_2,out}\right)x_{H_2}{C}_{H_2}\left(T_{fuel,out}-T_{fuel,in}\right)\\ & +\left(M_{H_{2}O,in}+M_{H_{2}O,out}\right)x_{H_{2}O}{C}_{H_{2}O}\left(T_{fuel,out}-T_{fuel,in}\right)\end{array};$$

  (28)
• Fuel:

  $$q_{fuel,net}=m_{fuel}{C}_{fuel}\frac{d{T}_{fuel}}{dt};$$

  (29)
• Air between the cell tube and the AST:

  $$q_{air,ann,net}=m_{air,ann}{C}_{air}\frac{d{T}_{air,ann}}{dt};$$

  (30)
• AST:

  $$q_{AST,net}=m_{AST}{C}_{AST}\frac{d{T}_{AST}}{dt},$$

  (31)

  $$q_{AST,conv,outer}=h_{AST,outer}{A}_{AST,outer}\left(T_{air,cell}-T_{AST}\right),$$

  (32)

  $$q_{AST,conv,inner}=h_{AST,inner}{A}_{AST,inner}\left(T_{AST}-T_{air,AST}\right),$$

  (33)

  $$q_{flow,air,AST}=m_{air}{C}_{air}\left(T_{air,AST,in}-T_{air,AST,out}\right);$$

  (34)
• Air in the AST:

  $$q_{air,AST,net}=m_{air,AST}{C}_{air}\frac{d{T}_{air,AST}}{dt}.$$

  (35)

Based on the formulation of (21)–(35), the energy balance model for the SOFC was developed.

2.2. Control Problems

The main propose of the controller design was to ensure that the output voltage V_out and the fuel utilization FU follow their reference values and maintain their stability as strongly as possible. The main problems in the controller design are listed as follows:

• High coupling: The SOFC is a coupled system, which makes its controller design challenging and complicated;
• Frequent disturbance: Due to the frequent change of the load, a continual disturbance exists in the system, which requires the controller to have strong robustness and accurate control performance;
• Conflicting objectives: The accuracy of the control performance of V_out and FU are the two conflicting objectives. However, to improve the power quality of the SOFC, both objectives need to be optimized at the same time.

3. Controller Design

3.1. Model Identification

Since the SOFC model is a multi-input–multi-output (MIMO) system, the load resistance R was chosen as the disturbance of the system, the fuel flow m_H2 and the cathode pressure P_c were chosen as the control variables, and the output voltage V_out and fuel utilization FU were chosen as the controlled variables. Thus, the linearized model can be identified as a 3 × 2 MIMO system:

$$\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right]=\left[\begin{array}{cc}{G}_{11}& G_{12}\\ G_{21}& G_{22}\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]+\left[\begin{array}{c}{G}_{\mathrm{d}1}\\ G_{\mathrm{d}2}\end{array}\right]d,$$

(36)

where y₁ is the output voltage V_out, y₂ is the fuel utilization FU, u₁ is the fuel flow m_H2, u₂ is the cathode pressure P_c, d is the load resistance R, G₁₁ represents the transfer function from m_H2 to V_out, G₁₂ represents the transfer function from m_H2 to FU, G₂₁ represents the transfer function from P_c to V_out, G₂₂ represents the transfer function from P_c to FU, G_d1 represents the transfer function from R to V_out, and G_d2 represents the transfer function from R to FU.

To identify these transfer functions, a step response of this system was analyzed. First, a step signal was set on m_H2, while P_c and R were kept as constants, and the step responses of V_out and FU were recorded. Second, a step signal was set on P_c, while m_H2 and R were kept as constants, and then the step responses of V_out and FU were recorded. Third, a step signal was set on R, while P_c and m_H2 were kept as constants, and the step responses of V_out and FU were recorded. Finally, based on the step response results, the transfer functions can be readily identified as follows:

$$\begin{array}{l}{G}_{11}=\frac{528.942s+2594.573}{s^2+0.627s+0.064}\\ G_{21}=\frac{0.108s+1.024\times {10}^{-4}}{s^2+0.118s+1.119\times {10}^{-4}}\\ G_{12}=\frac{-271.283s-0.286}{s^2+0.107s+1.116\times {10}^{-4}}\\ G_{22}=\frac{-8.998\times {10}^{-6}s-5.227\times {10}^{-5}}{s^2+0.209s+0.012}\\ G_d{}_1=\frac{0.591s+5.122\times {10}^{-4}}{s^2+0.113s+9.991\times {10}^{-5}}\\ G_d{}_2=\frac{-0.015s-1.306\times {10}^{-5}}{s^2+0.121s+1.050\times {10}^{-4}}\end{array}$$

(37)

The step response curves are drawn in Figure 2. As we can tell from the results, the step response curves of the identified model are the same as those of the original model, and thus the identified model is able to be used for control design.

3.2. Relative Gain Array (RGA) Paring

To analyze the relationship between the input and output variables of a system, we applied the relative gain array (RGA) defined by Skogestad and Postlethwaite in [27]. First, by setting the time s to 0, we got the steady-state matrix of G₁₁, G₁₂, G₂₁, and G₂₂:

$$A=\left[\begin{array}{cc}{G}_{11}{|}_{s\to 0}& G_{12}{|}_{s\to 0}\\ G_{21}|{}_{s\to 0}& G_{22}|{}_{s\to 0}\end{array}\right]=\left[\begin{array}{cc}4.0763\times {10}^4& 0.9148\\ -2.5607\times {10}^3& 0.0045\end{array}\right].$$

(38)

Then, the RGA matrix was calculated:

$$RGA=A{\left(A^{-1}\right)}^T=\left[\begin{array}{cc}{\lambda }_{11}& {\lambda }_{12}\\ {\lambda }_{21}& {\lambda }_{22}\end{array}\right]=\left[\begin{array}{cc}0.0729& 0.9271\\ 0.9271& 0.0729\end{array}\right].$$

(39)

The results show that the output voltage V_out is affected more by the cathode pressure P_c, while the fuel utilization FU is affected more by the fuel flow m_H2. Thus, V_out can be controlled by P_c easily, and FU can be controlled by m_H2 easily.

3.3. Controller Design

In this paper, two PID controllers were applied to control m_H2 and P_c independently. The formula for a PID controller is shown as

$$G\left(s\right)=k_p+k_i\frac{1}{s}+k_d\frac{T_d}{1+T_d\frac{1}{s}},$$

(40)

where k_p, k_i, and k_d are the proportional gain, integral gain, and derivative gain, respectively. In previous research, the main work has been to search for the optimal values of k_p, k_i, and k_d [17,19,28]. However, the filter coefficient N is also of great importance, and an appropriate filter coefficient can reduce the output oscillation and improve the control performance. As a result, the filter coefficient was also studied for this paper, and there were eight parameters that needed to be decided on, which were the parameters of PID1 (k_p1, k_i1, k_d1, N₁) and the parameters of PID2 (k_p2, k_i2, k_d2, N₂).

The structure of the control system is shown in Figure 3, and the parameters of the controllers were initialized by the pidTuner Toolbox in Matlab as

$$\{\begin{array}{l}{k}_{p1}=-1.646\times {10}^{-4}\\ k_{i1}=-1.360\times {10}^{-4}\\ k_{d1}=-5.848\times {10}^{-6}\\ T_{d1}=16.970\end{array},\{\begin{array}{l}{k}_{p2}=1.572\\ k_{i2}=0.326\\ k_{d2}=-0.067\\ T_{d2}=1.012\end{array}.$$

(41)

4. Solution Method

4.1. Introduction of Multiobjective Optimization based on a Genetic Algorithm

Multiobjective optimization usually involves two or more objective functions, and it is widely used in the field, where objectives are conflicting and need to be optimized at the same time. As is shown in Figure 4, the optimization result is a set of Pareto optimal solutions, where none of the objectives can be promoted without deteriorating other objectives. A multiobjective optimization problem (MOP) can be described as

$$\begin{array}{c}\min J\left(x\right)=[j_1\left(x\right),j_1\left(x\right),\cdots ,j_n\left(x\right)]\\ x\in X\end{array},$$

(42)

where J is the vector of objective functions, n is the number of objective functions, x represents the decision vector, and X is the constraint of the decision vector x. The objective functions are designed by the decision-maker to judge the quality of the solution, and in the field of industrial control, objective functions usually consist of an IAE, overshoot, and settling time.

To solve optimization problems, many methods have been developed, and the GA proposed by Holland in 1975 [29] is one of the most popular ones. Meanwhile, in 1994, Srinivas and Deb [30] developed a new algorithm based on a GA called a nondominated sorting genetic algorithm (NSGA), and then in 2000, an improved version of the NSGA (called the NSGA-II) was developed by Deb et al. [31]. With the merits of low computational complexities and the ability of parallel computation, it is widely used for MOPs [21,32]. Afterwards, in order to improve the distribution density of the Pareto front and avoid the problem of premature ripening, an improved NSGA-II was developed by introducing the parameter Pareto faction x, which represents the ratio of the optimal solutions to the total population (Equation (43)):

$$x=\frac{N_{\mathrm{Pareto}}}{N},$$

(43)

where N represents the population size, and N_Pareto represents the number of Pareto optimal solutions. A flow chart of this algorithm is shown in Figure 5.

In this section, the proposed improved NSGA-II is applied to find the optimal sets of control parameters for the SOFC.

4.2. Optimization of the Control Parameters of the SOFC

As was mentioned in Section 3, there were eight parameters in the two PID controllers that needed to be optimized: k_p1, k_i1, k_d1, N₁, k_p2, k_i2, k_d2, and N₂. Thus, the decision vector could be described as x = (k_p1, k_i1, k_d1, N₁, k_p2, k_i2, k_d2, N₂). The integral absolute error (IAE) of the output voltage and fuel utilization were employed as two objective functions, and the MOP could be formulated as

$$\min F=\left[\begin{array}{cc}IA{E}_1\left(x\right)& IA{E}_2\left(x\right)\end{array}\right],$$

(44)

where

$$IA{E}_1={\displaystyle \underset{0}{\overset{T}{\int }}\left(\left|V_{\mathrm{out},\mathrm{ref}}-V_{\mathrm{out}}\right|\right)dt},$$

(45)

$$IA{E}_2={\displaystyle \underset{0}{\overset{T}{\int }}\left(\left|F{U}_{\mathrm{ref}}-FU\right|\right)dt},$$

(46)

and the constraint for each parameter is shown as

$$\{\begin{array}{l}{k}_{p1}\in \left(-{10}^{-2},-{10}^{-7}\right)\\ k_{i1}\in \left(-{10}^{-2},-{10}^{-7}\right)\\ k_{d1}\in \left({10}^{-8},{10}^{-4}\right)\\ T_{d1}\in \left({10}^{-3},100\right)\end{array},\{\begin{array}{l}{k}_{p2}\in \left({10}^{-3},10\right)\\ k_{i2}\in \left({10}^{-3},10\right)\\ k_{d2}\in \left(-1,-{10}^{-4}\right)\\ T_{d2}\in \left({10}^{-3},100\right)\end{array}.$$

(47)

The population size N was set as 50, the maximum generation Max was set as 50, the Pareto fraction x was set as 0.35, the mutation function was a constraint-dependent type, and the crossover probability was 0.8. The stopping criteria were designed as follows:

• The number of generations is higher than the maximum number of iterations;
• The average relative change in the spread of the Pareto solutions over 100 generations is less than 10⁻⁴, and the spread is smaller than the average spread over the last 100 generations.

For each generation, four load disturbances occurred in each simulation process, and two objective functions, IAE₁ and IAE₂, of 50 individuals were calculated based on the simulation results.

The distribution of the final population in the objective function space is shown in Figure 6, which clearly illustrates the conflicting relationship between IAE₁ and IAE₂. Any decrease in IAE₂ led to an increase in IAE₁. In other words, any improvement in the control performance of fuel utilization led to a deterioration in the control performance of the output voltage. Five typical points were chosen from the Pareto front to get the optimal parameters of the controllers, and the value of each parameter and the objective function are displayed in

Table 1. The value of the control parameters and objective functions at the five optimal points.

Control Parameters	A	B	C	D	E
kp1	−0.00753	−0.00666	−0.00708	−0.00677	−0.00733
ki1	−0.00981	−0.00983	−0.00985	−0.00999	−0.00998
kd1	8.18 × 10−5	7.68 × 10−5	5.40 × 10−5	6.58 × 10−5	5.88 × 10−5
Td1	56.01119	55.96752	68.6119	53.9818	48.92057
kp2	9.950384	9.953341	7.995275	9.506123	8.360098
ki2	9.977686	9.975138	9.200353	2.342223	1.672613
kd2	−0.92397	−0.90781	−0.51078	−0.63197	−0.50791
Td2	69.30712	69.32299	51.57974	69.18918	63.97332
IAE1	1.203	1.205	1.481	3.081	4.313
IAE2	0.0134	0.0113	0.0111	0.011	0.0108

. In Table 1, point A possessed the lowest IAE₁ and the highest IAE₂, while point E possessed the highest IAE₁ and the lowest IAE₂. As we can tell, there were trade-offs between the two objectives, and to decide the final optimal solution, a decision-making process should be executed based on the decision-maker’s insight.

Figure 7 denotes the variation of the objective functions with the generation. As we can tell from the result, IAE₁ decreased from 22.75 to 1.20 and converged at about 1.197 after 30 generations, while IAE₂ decreased from 0.8 to 0.011 and converged at about 0.0108 after 25 generations. Both IAE₁ and IAE₂ were greatly minimized after about 30 generations, and thus the optimization was successful.

5. Results and Discussion

In this section, to analyze the control performance of the five sets of optimal control parameters and the initial control parameters, a simulation was performed based on the linearized model implemented, and the results were compared and discussed.

5.1. Comparison of the Optimal Points

In this part, a step signal was added to the load resistance, and the disturbance responses of the controllers on the five optimal points and the original controller were compared. The best optimal points were chosen based on the comparison results.

5.1.1. Simulation Results

During the simulation process, a step signal was added to the load resistance to make it change from 5 Ω to 6 Ω at 20 s, while the reference values of V_out and FU were 168 V and 80%, respectively. The disturbance responses of the output power P_out, the output voltage, the fuel utilization, the fuel flow, and the cathode pressure are shown in Figure 8, and detailed data from the simulation are shown in

Table 2. Detailed data from the simulation results.

Output Variables	Index	Before Optimization	Optimal Point A	Optimal Point B	Optimal Point C	Optimal Point D	Optimal Point E
Output Voltage	Overshoot (%)	0.0227	0.0113	0.0111	0.0211	0	0
	Settling time (s)	38.231	4.293	4.356	4.903	10.865	17.108
	IAE1	11.7506	0.2897	0.2924	0.3495	0.8972	2.0206
Fuel Utilization	Overshoot (%)	0	0	0	0	0	0
	Settling time (s)	4.832	1.225	0.927	0.0937	0.0652	0.0623
	IAE2	0.2611	0.0037	0.0028	0.0025	0.0018	0.0017

.

5.1.2. Discussion

From the simulation results, we can observe that compared to the initial controller before optimization, the controllers tuned to the five optimal points all had better control performance.

In Figure 8b, the controllers on the five optimal points all had faster dynamic responses on P_out compared to the initial controller. At optimal point A, the speed of the dynamic response was the fastest, while optimal point E possessed the slowest dynamic response. In Figure 8c,d, the controllers on the five optimal points all had faster control effects on P_c and m_H2.

In Figure 8e, it is clear that the control performance of V_out could be dramatically improved by applying the controllers on the five optimal points, and the settling time and IAE of the output voltage gradually increased from optimal point A to E, while the optimal point C possessed the biggest overshoot among all the optimal points (0.0211%). In Figure 8f, both the settling time and IAE of the response of FU on the five optimal points were better than those before optimization. Among all of the optimal points, the controller on point A possessed the longest settling time (1.225 s), while the settling time for the controller before optimization was 4.832 s. In addition, the controller on optimal point E possessed the best control performance with the shortest settling time, 0.062 3 s, and the lowest IAE, 0.0017.

In conclusion, the controllers on the five optimal points all had a smaller overshoot, lower settling times, and a lower IAE than the initial controller. Based on the simulation results, the control parameters on optimal point D were chosen as the best solution in order to reduce the control action and improve the control performance of both the output voltage and fuel utilization.

Additionally, there were some defects in the optimization process. Compared to the improved NSGA-II we applied in this paper, the NSGA-III proposed by Deb and Jain [33] in 2014 is more proper for an MOP with a Pareto front, whose objective values may have a different scale. Moreover, the population size and execution times applied were the default values in the Matlab Optimization Toolbox, which could be optimized to save on optimization time.

5.2. Continuous Disturbance Response Simulation

In this part, continuous step signals were added to the load resistance and the current to simulate the realistic working conditions of an SOFC. The disturbance responses of the optimal controller based on optimal point D and the original controller were compared, and conclusions were drawn based on the simulation results.

5.2.1. Resistance Disturbance Response Simulation

During the simulation process, a continuous step signal was added to the load resistance to make it change from 5 Ω to 5.5 Ω, 6 Ω, 6.5 Ω, and 5 Ω at 200 s, 400 s, 600 s, and 800 s, respectively. The reference values of V_out and FU were kept at 168 V and 80%, respectively. The disturbance responses of the output power P_out, the output voltage, the fuel utilization, the fuel flow, and the cathode pressure are shown in Figure 9, and detailed data from the simulation are shown in

Table 3. Detailed data from the simulation results.

Output Variables	Number	Optimal Controller		Original Controller
		Overshoot (%)	Settling Time (s)	Overshoot (%)	Settling Time (s)
Output Voltage	1	0	9.6	0.012	30.3
	2	0	9.3	0.007	24.5
	3	0	8.5	0	18.1
	4	0	16.4	0.060	46.2
Fuel Utilization	1	0	1.4	0	5.4
	2	0	1.5	0	5.2
	3	0	1.1	0	4.8
	4	0	2.4	0	6.8

.

5.2.2. Current Disturbance Response Simulation

During the simulation process, a continuous step signal was added to the current to make it change from 33 A to 27 A, 29 A, 31 A, and 33 A at 200 s, 400 s, 600 s, and 800 s, respectively. The reference values of V_out and FU were kept at 168 V and 80%, respectively. The disturbance responses of the output power P_out, the output voltage, the fuel utilization, the fuel flow, and the cathode pressure are shown in Figure 10, and detailed data from the simulation are shown in

Table 4. Detailed data from the simulation results.

Output Variables	Number	Optimal Controller		Original Controller
		Overshoot (%)	Settling Time (s)	Overshoot (%)	Settling Time (s)
Output Voltage	1	0	17.3	0.018	37.6
	2	0	11.9	0.009	38.8
	3	0	12.7	0.012	40.5
	4	0	12.1	0.015	50.3
Fuel utilization	1	0	2.2	0	6.5
	2	0	1.5	0	5.6
	3	0	1.4	0	5.2
	4	0	1.8	0	6.7

.

5.2.3. Discussion

From the simulation results, it is evident that compared to the original controller, the optimal controller based on the controller parameters on optimal point D had better control performance.

In Figure 9b and Figure 10b, the step responses of the optimal controller possessed faster dynamic responses than the original one, which led to fewer fluctuations in the output power. Figure 9c,d and Figure 10c,d, also illustrate that the optimal controller had more rapid dynamic responses to cathode pressure and fuel flow, which permitted faster control action in the controller.

From Figure 9e and Figure 10e, it is clear that the control performance of the optimal controller was much better than the original controller from the perspective of the overshoot and settling time. The overshoot of all of the disturbance responses of the optimal controller was 0, whereas the original controller had a slight overshoot. As for the settling time, the settling time of the optimal controller under resistance disturbance varied from 8.5 s to 16.4 s, while that of the original controller under resistance disturbance varied from 18.1 s to 46.2 s. In addition, the settling time of the optimal controller under current disturbance varied from 12.1 s to 17.3 s, while that of the original controller under resistance disturbance varied from 37.6 s to 50.3 s.

In Figure 9f and Figure 10f, although the overshoot of both the original controller and optimal controller was 0, the settling time of the optimal controller was much shorter due to the rapid control action. The settling time of the optimal controller under resistance disturbance varied from 1.1 s to 2.4 s, and as for the original controller, it varied from 4.8 s to 6.8 s. The settling time of the optimal controller under current disturbance varied from 1.8 s to 2.2 s, and the original controller varied from 5.2 s to 6.7 s.

In conclusion, from the simulation results of the continuous disturbance response, it is clear that the optimal controller on optimal point D had the merits of stability, a smaller overshoot, and reduced settling time, and consequently, it is the ideal controller for the SOFC model under consideration.

6. Conclusions

An SOFC is an efficient energy conversion device that has been widely used. To deal with the characteristics of couplings of multivariables, nonlinearity of the model, and frequent load disturbances, this paper employed a GA-based optimization algorithm to find the optimal parameters for the PID controllers of an SOFC. In this paper, a 3 × 2 MIMO system was identified based on the SOFC model, and two PID controllers were designed to separately control V_out and FU. Then, the optimization objectives were designed, and a GA was applied for multiobjective optimizations of the controller parameters for the two controllers. Finally, a set of optimal points on the Pareto front were found, and the corresponding control parameters of the five points were applied to control the linearized SOFC model. The results showed that the controller parameters at the five optimal points resulted in better control performance compared to before optimization in terms of the IAE, the overshoot, and the settling time. The parameters on optimal point D were considered to be the best for their smooth control action. Therefore, PID controllers with the controller parameters on optimal point D are an ideal choice for SOFC control.