Adaptive Tracking Constrained Controller Design for Solid Oxide Fuel Cells Based on a Wiener-Type Neural Network

Abstract

In order to solve the control problem of the solid oxide fuel cell(SOFC), a novel adaptive tracking constrained control strategy based on a Wiener-type neural network is proposed in this paper. The working principle of SOFC is introduced, and the dynamical model of SOFC is studied. Besides, a Wiener model formulation for SOFC is proposed to approximate the nonlinear dynamics of the system, and an adaptive Wiener model identification method is utilized to identify the parameters of the model. Moreover, an adaptive exponential PID controller is designed to keep the stack output voltage stable. Meanwhile, the saturation problem is considered in the paper including input magnitude and rate constraints. Additionally, an anti-windup compensator is employed to eliminate the abominable influence of the saturation problem. Then, the stability of the control plant is analyzed and proven via the Lyapunov function. Finally, the simulation based on the MATLAB/Simulink environment is carried out, and the conventional PID controller is added and simulated as a contrast to verify the control performance of the proposed control algorithm. The results indicate that the proposed control algorithm possesses favorable control performance when dealing with nonlinear systems with complex dynamics.

1. Introduction

In recent years, with the gradual scarcity of non-renewable energy resources, people have been eager to find alternatives to fossil fuel. Against this background, many researchers have considered fuel cell electricity generation with higher energy efficiency, which produces water, heat, less noise, no radiation and negligible amounts of pollutants [1,2,3,4,5]. The solid oxide fuel cell (SOFC) is a type of fuel cell that uses solid oxide material as the electrolyte to send negative oxygen ions from the cathode to the anode, where a chemical reaction related to the electrochemical oxidation occurs [6,7,8,9,10]. SOFCs have many advantages, such as flexibility in selecting fuel, long-term stability and no need for an expensive catalyst. However, some serious problems due to sluggish dynamics, nonlinearity, fuel starvation, external disturbance, unavailability of some key measurements and tight operating constraints often occur in the actual control process of an SOFC [11]. Therefore, the control strategy of the SOFC system is still a challenging problem.

The complex dynamical system of an SOFC is often subject to imprecise models and a priori unknown disturbances. Therefore, some model-based control methods cannot satisfy the control objective due to the model-reality discrepancy [12,13]. To obtain an approximate model of the SOFC under unknown parameters, many algorithms based on the neural network model [14,15,16] have been proposed. In this paper, a Wiener-type neural network is employed to achieve the identification of the model of SOFC. On the one hand, neural networks are widely used in identification and control of nonlinear systems with complex dynamics due to their ability to approximate any nonlinear function. They have several advantages, for example, the unique learning ability of a neural network makes it able to adapt to the changes of the system and environment. The powerful parallel computing ability makes it possible to carry out a large number of complex operations quickly. Besides, the multi-input multi-output model structure makes it possible to facilitate identification and control of multivariable systems. In [14], a fuzzy neural network controller was proposed for the control of active and reactive power of an SOFC-based system. In [15], an artificial neural network was used to establish the model of SOFC and predict the performance and parameters of SOFC. In [16], a nonlinear predictive controller based on an improved radial basis function neural network identification model was designed to guarantee the fuel utilization to operate within a safe range. On the other hand, the Wiener model composed of a linear dynamic system and a static nonlinear system in series has a clear internal structure and can be used to express the nonlinear characteristics of common systems in most cases. Besides, the model is often employed based on some intelligence algorithms like a neural network [17,18]. It has attracted wide attention for a long time and has been applied in many research fields. In [17], the authors gave a method of the Wiener identification model in series connected by a linear dynamic neuron and static BP neural network. In [18], a nonlinear Wiener model combined with the Laguerre function and a BP network was presented, which had strong robustness due to the information of the model being non-essential.

In addition, the problem of input constraints, which can impact the system performance and stability, is ignored by many researchers. Notice that in terms of preventing fuel from being overused in an SOFC system, the expected scope of fuel utilization must be maintained from 0.7–0.9, which should further take the input constraints into consideration. Then, various constraint problems have been considered in closed-loop control systems [19,20,21,22,23]. At present, most of the constrained control approaches for the SOFC system focus on model predictive control (MPC) [5,24,25,26,27]. In [26], the researchers used a new offset-free input to design a fuzzy controller. A predictive control method based on the fuzzy Hammerstein model was shown in [27]. However, MPC has some typical drawbacks in actual industrial processes. Firstly, it is difficult to obtain a reliable first-principles model for an object with highly nonlinear characteristics. Secondly, the stability analysis of MPC is often hard to obtain and yet can be easily achieved via other control methods like model-free adaptive control [22,23]. Finally, the cost of calculation for MPC is quite large [26,28,29]. To avoid the disadvantages of MPC, an adaptive tracking constrained control method is proposed in this paper. The proposed control method depends on the modified PID control algorithm, whose exponential parameters can be estimated by an adaptive law based on a back-propagation training algorithm. Moreover, the complex dynamics of SOFC can be approximated via the Wiener-type neural network, and stability analysis of the proposed control algorithm can be easily completed by the Lyapunov function. Furthermore, its calculation is reduced in contrast with MPC.

In this paper, we first present an innovative, novel adaptive constrained control design based on the Wiener-type neural network identification for SOFCs. There are various amplitudes and rate constraints of the control inputs in the SOFC control system; based on the online Wiener-type neural network model identification and exponential PID control, an innovative dynamic anti-windup compensator is added to eliminate the negative effects of control input constraints. Furthermore, a stability analysis of the SOFC is given for the presented adaptive constrained approach via the Lyapunov function. Simulations of the SOFC prove that our adaptive tracking constrained control strategy based on Wiener model identification is highly effective. The main innovation points in this paper are as follows:

• An innovative method based on the Wiener-type neural network model can be used to identify the system dynamics of an SOFC.
• A new dynamic anti-windup compensator combined with an exponential PID control method is proposed to eliminate the influence of amplitude and rate constraints for control input in the SOFC system.
• The Lyapunov function is used to verify that the SOFC closed-loop control system is stable.

The rest of this paper is as follows. The problem formulation of the SOFC voltage control system is illustrated in Section 2. In Section 3, the description of the constrained adaptive tracking control based on the Wiener model is presented. Simulation results for the SOFC voltage control system are given in Section 4. Some conclusions for the whole paper are drawn in Section 5.

2. Problem Formulation For SOFC

According to the previous research [7,8,9,10], the working principles of fuel cells are essentially similar. In short, the SOFC system contains the fuel process part and the fuel cell stack. Functions of the SOFC are the same as those of the inverse device for water electrolysis. The general structure of the SOFC is depicted in Figure 1. The single cell for the SOFC includes several parts, such as the anode, the cathode and the solid oxide electrolyte. Fuel oxidation and oxidant reduction occur in the anode and cathode, respectively. At the same time, the catalyst, which can accelerate the electrochemical reaction of the electrode, is widely available in the anode and cathode. Based on the above introduction of the SOFC, it can be speculated that the actual effectiveness of the SOFC is equivalent to a DC voltage source in that the SOFC system can provide electricity continuously via the electrochemical reactions. Comparing the structures of the SOFC and DC source, it is indicated that the anode of the SOFC is equivalent to the negative power supply of the DC source, while the cathode of the SOFC is equivalent to the positive power supply of the DC source.

The basic process of the chemical reaction in an SOFC is introduced in Figure 1. The fuel gas, similar to hydrogen (${\mathrm{H}}_2$), methane (${\mathrm{CH}}_4$) and city gas (CO), must be supplied on the anode side continuously. Then, the surface of the anode, which has a porous structure, adsorbs the fuel gas and sends it to the insideof the anode and electrolyte. ${\mathrm{O}}_2$ can be transformed into ${\mathrm{O}}^{2-}$ when the catalytic action inside the cathode gives ${\mathrm{O}}_2$ an electron. Due to the impact of the chemical potential, ${\mathrm{O}}^{2-}$ combines with the solid oxygen ion conductor, which can be realized as an electrolyte. Under the influence of the concentration gradient, ${\mathrm{O}}^{2-}$ can react with the fuel gas when ${\mathrm{O}}^{2-}$ diffuses into the inside of the electrolyte and anode. Then, all lost electrons return to the cathode to combine with ${\mathrm{O}}_2$. Overall, this is a cyclic process.

The main electrochemical processes of the SOFC can be described as follows:

$$\begin{array}{cc}\hfill \mathrm{Anode}:& {\mathrm{H}}_2+{\mathrm{O}}^{2-}\to {\mathrm{H}}_2\mathrm{O}+2{\mathrm{e}}^{-}\hfill \\ \hfill \mathrm{Cathode}:& {\mathrm{O}}_2+4{\mathrm{e}}^{-}\to 2{\mathrm{O}}^{2-}\hfill \end{array}$$

In this paper, we adopt the widely-accepted nonlinear dynamical model of the SOFC, which was presented in [27,30]. This SOFC model is depicted in Figure 2. Taking advantage of Nernst’s equation with the consideration of ohmic, concentration and activation losses (i.e., ${\eta }_{ohmic}$, ${\eta }_{conc}$ and ${\eta }_{act}$), the function of the stack output voltage $V_{dc}$ associated with the natural gas flow rate $q_f$ can be defined as:

$$\begin{array}{c}\hfill V_{dc}=V_0-{\eta }_{ohmic}-{\eta }_{conc}-{\eta }_{act}\end{array}$$

(1)

therein

$$\begin{array}{c}\hfill V_0=N_0\left[E_0+\frac{R_0{T}_0}{2F_0}\ln \frac{p_{{\mathrm{H}}_2}\sqrt{p_{{\mathrm{O}}_2}/101,325}}{p_{{\mathrm{H}}_2\mathrm{O}}}\right]\end{array}$$

(2)

$$\begin{array}{c}\hfill p_{{\mathrm{H}}_2}=\frac{1}{K_{{\mathrm{H}}_2}(1+{\tau }_{{\mathrm{H}}_2}s)}\left(\frac{1}{1+{\tau }_{f}s}{q}_f-2K_{r}I\right)\end{array}$$

(3)

$$\begin{array}{c}\hfill p_{{\mathrm{O}}_2}=\frac{1}{K_{{\mathrm{O}}_2}(1+{\tau }_{{\mathrm{O}}_2}s)}\left(\frac{1/{\tau }_{\mathrm{H}-\mathrm{O}}}{(1+{\tau }_{f}s)}{q}_f-K_{r}I\right)\end{array}$$

(4)

$$\begin{array}{c}\hfill p_{{\mathrm{H}}_2\mathrm{O}}=\frac{2}{K_{{\mathrm{H}}_2\mathrm{O}}(1+{\tau }_{{\mathrm{H}}_2\mathrm{O}}s)}{K}_{r}I\end{array}$$

(5)

$$\begin{array}{c}{\eta }_{ohmic}=Ir,{\eta }_{conc}=\partial +\beta \ln I,\hfill \\ {\eta }_{act}=-\frac{R_0{T}_0}{2F_0}\ln I\left(1-\frac{I}{I_L}\right)\hfill \end{array}$$

(6)

For the nonlinear dynamic system of the SOFC, the parameters and corresponding values are shown in

Table 1. Parameters of the SOFC system model.

Parameter	Value	Parameter	Value
$T_0$	1273	${\tau }_{{\mathrm{H}}_2}$	26.1
$F_0$	96,485	${\tau }_{{\mathrm{H}}_2\mathrm{O}}$	78.3
$R_0$	8.314	${\tau }_{{\mathrm{O}}_2}$	2.91
$E_0$	1.18	${\tau }_{\mathrm{H}-\mathrm{O}}$	1.145
$N_0$	384	r	0.126
$K_r$	$0.996\times {10}^{-3}$	${\tau }_f$	5
$K_{{\mathrm{H}}_2}$	$8.32\times {10}^{-6}$	∂	0.05
$K_{{\mathrm{H}}_2\mathrm{O}}$	$2.77\times {10}^{-6}$	$\beta$	0.11
$K_{{\mathrm{O}}_2}$	$2.49\times {10}^{-5}$	$I_L$	800

[24].

In the design process of this dynamic control system for the SOFC, we mainly consider three aspects, which are described as follows:

• It is found that the current load I is related to the stack output voltage $V_{dc}$ based on the above dynamic model (see Figure 2), and changes of the current load I can influence the output voltage $V_{dc}$ directly. Therefore, one of the main control objectives is to eliminate the negative impact due to changes of the current load to obtain a stable output voltage.
• Based on the stationary voltage-current characteristics of the SOFC (see Figure 3), it is shown that the behavior of the SOFC reflects the typical nonlinear characteristic over a wide operating regime. In addition, the nonlinear characteristic of system responses in the SOFC system at low- and high-current loads is obvious. Even under the condition of an overloaded current, the operating stack voltage may deteriorate rapidly, which can also reflect the strong nonlinear characteristic of SOFCs. Thus, the nonlinear control strategy must be better than the linear method for SOFCs. In the design of the nonlinear controller, $q_f$, $V_{dc}$ and I need to be defined as the control input u, control output y and measurable disturbance d, respectively. This input-output relationship can be redefined as an unknown discrete-time nonlinear autoregressive exogenous model (NARX):

  $$\begin{array}{c}y\left(k\right)=f\left(y(k-1),\cdots ,y(k-n_a),u\left(k\right),u(k-1),\cdots ,\right.\\ \left.u(k-n_b),d(k-1),\cdots ,d(k-n_c)\right)\end{array}$$

  (7)

  where $n_a$, $n_b$ and $n_c$ denote the orders of dynamics. Moreover, $n_a\ge n_b$, $n_a\ge n_c$.
• In practice, the control input $q_f=u$ does not change rapidly in a short time with inertia existing in the SOFC control system. For this SOFC system, the control inputs must be constrained by the amplitude and rate conditions, which can be written as:

  $$\begin{array}{c}{q}_{f\min }\le u\left(t\right)\le q_{f\max },{\dot{q}}_{f\min }\le \dot{u}\left(t\right)\le {\dot{q}}_{f\max }\hfill \end{array}$$

  (8)
  In addition to the constraints of the control input, fuel utilization is another important operating condition that impacts the performance of an SOFC. Thus, the fuel utilization needs to be maintained within a limited scope, which can be described as:

  $$\begin{array}{c}\hfill \rho =\frac{q_{{\mathrm{H}}_2}^{in}-q_{{\mathrm{H}}_2}^o}{q_{{\mathrm{H}}_2}^{in}}=\frac{q_{{\mathrm{H}}_2}^r}{q_{{\mathrm{H}}_2}^{in}}=\frac{2K_{r}I}{q_{{\mathrm{H}}_2}^{in}}\end{array}$$

  (9)
  The expected utilization must be from ${\rho }_{\min }=0.7$ to ${\rho }_{\max }=0.9$. The hydrogen flow should be within $\left[q_{{\mathrm{H}}_2\min }^{in},q_{{\mathrm{H}}_2\max }^{in}\right]$ from Equation (9) to maintain safe utilization. Then, the result $\left[q_{{\mathrm{H}}_2\min }^{in},q_{{\mathrm{H}}_2\max }^{in}\right]$ can be employed to compute the constrained condition $\left[{\acute{q}}_{f\min },{\acute{q}}_{f\max }\right]$ of the fuel flow. Finally, the ultimate dynamic amplitude constraint $\left[{\overline{q}}_{f\min },{\overline{q}}_{f\max }\right]$ can be obtained as:

  $$\begin{array}{cc}\hfill {\overline{q}}_{f\min }& =\max \left(q_{f\min },{\acute{q}}_{f\min }\right)\hfill \\ \hfill {\overline{q}}_{f\max }& =\min \left(q_{f\max },{\acute{q}}_{f\max }\right)\hfill \end{array}$$

  (10)
  Therefore, the following relationship can be achieved as:

  $$\begin{array}{c}{\overline{q}}_{f\min }\le u\left(k\right)\le {\overline{q}}_{f\max }\hfill \end{array}$$

  (11)
  Due to the final dynamic magnitude constraint, the change rate of $u\left(t\right)$ can be approximated via the first-order Euler approximation, which is written as:

  $$\dot{u}\left(t\right)=\frac{u\left(k\right)-u(k-1)}{T_s}$$

  (12)

  where $T_s$ denotes the sampling time. According to (8) and (12), the following function can be derived as:

  $$T_s{\dot{q}}_{f\min }+u(k-1)\le u\left(k\right)\le T_s{\dot{q}}_{f\max }+u(k-1)$$

  (13)
  From (11) and (13), the final unified input constraint for the SOFC is summarized as:

  $$\begin{array}{c}\hfill \max \left\{{\overline{q}}_{f\min },\left(T_s{\dot{q}}_{f\min }+u(k-1)\right)\right\}\le u\left(k\right)\le \min \left\{{\overline{q}}_{f\max },\left(T_s{\dot{q}}_{f\max }+u(k-1)\right)\right\}\end{array}$$

  (14)

Remark 1.

It is generally believed that optimal control is the unique comprehensive strategy to address the problem of multiple constraints for the SOFC due to the inherent capability of this algorithm. Note that MPC gives us a convenient way to solve the problems of various constraints, including amplitude, rate and higher-order saturations [31]. However, there is still a problem that cannot be neglected, i.e., the stability analysis of general optimal control for a nonlinear system with multiple constraints is difficult, and many optimal control methods require numerous computations, yielding an unsatisfactory performance.

3. Constrained Adaptive Tracking Control

3.1. Wiener Model Formulation for SOFC

The Wiener model can approach the nonlinear function arbitrarily, so it is designed to identify the dynamic model. In some literature, the general Wiener model can be composed of the linear dynamic part and the nonlinear static part [32]; this structure chart of the Wiener model is depicted in Figure 4.

For the SOFC control system (7), the linear dynamic block of the Wiener model can be defined as:

$$\begin{array}{c}x\left(k\right)=\frac{B\left(q^{-1}\right)}{A\left(q^{-1}\right)}u\left(k\right)+\frac{C\left(q^{-1}\right)}{A\left(q^{-1}\right)}d\left(k\right)\hfill \end{array}$$

(15)

where:

$$\begin{array}{c}A\left(q^{-1}\right)=1+a_1{q}^{-1}+a_2{q}^{-2}+\cdots +a_{n_a}{q}^{-n_a}\hfill \\ B\left(q^{-1}\right)=b_0+b_1{q}^{-1}+b_2{q}^{-2}+\cdots +b_{n_b}{q}^{-n_b}\hfill \\ C\left(q^{-1}\right)=c_1{q}^{-1}+c_2{q}^{-2}+\cdots +c_{n_c}{q}^{-n_c}\hfill \end{array}$$

Therein, $a_i$, $b_j$ and $c_l$ denote the parameters of the orders, $i=1,\cdots ,n_a$, $j=0,1,\cdots ,n_b$ and $l=1,\cdots ,n_c$.

In this paper, the polynomial forms are utilized to describe the dynamic characteristics of the nonlinear static block [33]. The polynomial function is defined as:

$$\begin{array}{c}\hfill y\left(k\right)=\sum _{m=1}^p{\sigma }_m{x}^m\left(k\right)\end{array}$$

(16)

where ${\sigma }_m$ is the parameter of the nonlinear static block and $m=1,\cdots ,p$; therein, m and p are the power of the intermediate variable $x\left(k\right)$ and the degree of the polynomial function, respectively.

3.2. Adaptive Wiener Model Identification

The scheme of the training Wiener model can easily identify the parameters $a_i$, $b_j$, $c_l$ and ${\sigma }_m$. The Wiener model output $\widehat{y}\left(k\right)$ is defined in the following form:

$$\begin{array}{c}\hfill \widehat{y}\left(k\right)=\sum _{m=1}^p{\widehat{\sigma }}_m{\widehat{x}}^m\left(k\right)\end{array}$$

(17)

In the Wiener model, the function of the hidden layer output is described as:

$$\begin{array}{ccc}\hfill \widehat{x}\left(k\right)& \hfill =& -{\widehat{a}}_1\widehat{x}(k-1)-{\widehat{a}}_2\widehat{x}(k-2)-\cdots -{\widehat{a}}_{n_a}\widehat{x}(k-n_a)\hfill \\ & & +{\widehat{b}}_{0}u\left(k\right)+{\widehat{b}}_{1}u(k-1)+\cdots +{\widehat{b}}_{n_b}u(k-n_b)\hfill \\ & & +{\widehat{c}}_{1}d(k-1)+{\widehat{c}}_{2}d(k-2)+\cdots +{\widehat{c}}_{n_c}d(k-n_c)\hfill \\ & \hfill =& -\sum _{i=1}^{n_a}{\widehat{a}}_i\widehat{x}(k-i)+\sum _{j=0}^{n_b}{\widehat{b}}_{j}u(k-j)+\sum _{l=1}^{n_c}{\widehat{c}}_{l}d(k-l)\hfill \end{array}$$

(18)

where the weights ${\widehat{a}}_i$, ${\widehat{b}}_j$, ${\widehat{c}}_l$ and ${\widehat{\sigma }}_m$ are related to the parameters $a_i$, $b_j$, $c_l$ and ${\sigma }_m$ in Equation (15) and (16), respectively.

Here, the negative gradient method can be employed to adjust the weights of the identified Wiener model. The identification error can be defined as:

$$\begin{array}{c}\hfill J_{W_1}=\frac{1}{2}{(y\left(k\right)-\widehat{y}\left(k\right))}^2=\frac{1}{2}{\widehat{e}}^2\left(k\right)\end{array}$$

(19)

Therein, $\widehat{e}\left(k\right)=y\left(k\right)-\widehat{y}\left(k\right)$. $W_1$ represents a set of ${\widehat{a}}_i$, ${\widehat{b}}_j$ ${\widehat{c}}_l$ and ${\widehat{\sigma }}_m$. The task obviously becomes an optimization problem, which makes the function Equation (19) about the weighting parameter $W_1$ minimum, which can be:

$$\begin{array}{c}\hfill \frac{\partial J_{W_1}}{\partial W_1}=-\widehat{e}\left(k\right)\frac{\partial \widehat{y}\left(k\right)}{\partial W_1}\end{array}$$

(20)

The general update rule can be defined in the following form:

$$\begin{array}{cc}\hfill W_1(k+1)& =W_1\left(k\right)+\Delta W_1\left(k\right)\hfill \\ & =W_1\left(k\right)+{\eta }_{w1}\widehat{e}\left(k\right)\frac{\partial \widehat{y}\left(k\right)}{\partial W_1}\hfill \end{array}$$

(21)

Therein, ${\eta }_{w1}$ represents the training rate of the weighting parameter $W_1$ in the Wiener model.

From Equations (17)–(19), the partial derivatives of output $\widehat{y}\left(k\right)$ to ${\widehat{\sigma }}_m$ and $\widehat{x}\left(k\right)$ in the Wiener model are written in the following form:

$$\begin{array}{c}\frac{\partial \widehat{y}\left(k\right)}{\partial {\widehat{\sigma }}_m}={\widehat{x}}^m\left(k\right)\hfill \\ \frac{\partial \widehat{y}\left(k\right)}{\partial \widehat{x}\left(k\right)}=\sum _{m=1}^{p}m{\widehat{\sigma }}_m{\widehat{x}}^{m-1}\left(k\right)\hfill \end{array}$$

(22)

Therein, $\widehat{x}\left(k\right)$ is the function of ${\widehat{a}}_i$, ${\widehat{b}}_j$ and ${\widehat{c}}_l$, and the partial derivatives of $\widehat{x}\left(k\right)$ about parameters ${\widehat{a}}_i$, ${\widehat{b}}_j$ and ${\widehat{c}}_l$ are described as:

$$\frac{\partial \widehat{x}\left(k\right)}{\partial {\widehat{a}}_i}=-\widehat{x}(k-i),\frac{\partial \widehat{x}\left(k\right)}{\partial {\widehat{b}}_j}=u(k-j),\frac{\partial \widehat{x}\left(k\right)}{\partial {\widehat{c}}_l}=d(k-l)$$

(23)

According to Equations (22) and (23), the partial derivatives of $\widehat{y}\left(k\right)$ about parameters ${\widehat{a}}_i$, ${\widehat{b}}_j$ and ${\widehat{c}}_o$ are defined in the following form:

$$\begin{array}{cc}\hfill \frac{\partial \widehat{y}\left(k\right)}{\partial {\widehat{a}}_i}& =\frac{\partial \widehat{y}\left(k\right)}{\partial \widehat{x}\left(k\right)}·\frac{\partial \widehat{x}\left(k\right)}{\partial {\widehat{a}}_i}=-\left(\sum _{m=1}^{p}m{\widehat{\sigma }}_m{\widehat{x}}^{m-1}\left(k\right)\right)\widehat{x}(k-i)\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill \frac{\partial \widehat{y}\left(k\right)}{\partial {\widehat{b}}_j}& =\frac{\partial \widehat{y}\left(k\right)}{\partial \widehat{x}\left(k\right)}·\frac{\partial \widehat{x}\left(k\right)}{\partial {\widehat{b}}_j}=\left(\sum _{m=1}^{p}m{\widehat{\sigma }}_m{\widehat{x}}^{m-1}\left(k\right)\right)u(k-j)\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill \frac{\partial \widehat{y}\left(k\right)}{\partial {\widehat{c}}_l}& =\frac{\partial \widehat{y}\left(k\right)}{\partial \widehat{x}\left(k\right)}·\frac{\partial \widehat{x}\left(k\right)}{\partial {\widehat{c}}_l}=\left(\sum _{m=1}^{p}m{\widehat{\sigma }}_m{\widehat{x}}^{m-1}\left(k\right)\right)d(k-l)\hfill \end{array}$$

(26)

Because of Equation (21), this update law about parameters ${\widehat{a}}_i$, ${\widehat{b}}_j$, ${\widehat{c}}_l$ and ${\widehat{\sigma }}_m$ can be rewritten as follows:

$$\begin{array}{cc}\hfill {\widehat{\sigma }}_m(k+1)& ={\widehat{\sigma }}_m\left(k\right)+{\eta }_{w1}\widehat{e}\left(k\right)\frac{\partial \widehat{y}\left(k\right)}{\partial {\widehat{\sigma }}_m}\hfill \\ \hfill {\widehat{a}}_i(k+1)& ={\widehat{a}}_i\left(k\right)+{\eta }_{w1}\widehat{e}\left(k\right)\frac{\partial \widehat{y}\left(k\right)}{\partial {\widehat{a}}_i}\hfill \\ \hfill {\widehat{b}}_j(k+1)& ={\widehat{b}}_j\left(k\right)+{\eta }_{w1}\widehat{e}\left(k\right)\frac{\partial \widehat{y}\left(k\right)}{\partial {\widehat{b}}_j}\hfill \\ \hfill {\widehat{c}}_l(k+1)& ={\widehat{c}}_l\left(k\right)+{\eta }_{w1}\widehat{e}\left(k\right)\frac{\partial \widehat{y}\left(k\right)}{\partial {\widehat{c}}_l}\hfill \end{array}$$

(27)

According to these update rules, the output $\widehat{y}\left(k\right)$ is determined. At the same time, the parameters ${\widehat{a}}_i$, ${\widehat{b}}_j$, ${\widehat{c}}_l$ and ${\widehat{\sigma }}_m$ in the Wiener model can be tuned. Based on the above research, the convergence criterion is proposed to select an adequate learning rate.

Theorem 1.

If ${\eta }_{w1}$ is chosen as:

$$\begin{array}{c}\hfill 0<{\eta }_{w1}<\frac{2}{{\overline{\mathrm{\Psi }}}_1^2}\end{array}$$

(28)

where ${\overline{\mathrm{\Psi }}}_1={\max }_k∥\frac{\partial \widehat{y}\left(k\right)}{\partial W_1}∥$, the modeling error $\widehat{e}\left(k\right)$ converges asymptotically to zero under the update law (27).

Proof of Theorem 1.

Choosing a Lyapunov function as $V_1\left(k\right)=\frac{1}{2}{\widehat{e}}^2\left(k\right)$, from (21), we get:

$$\begin{array}{cc}\hfill \Delta V_1\left(k\right)& =V_1(k+1)-V_1\left(k\right)\hfill \\ & =\Delta \widehat{e}\left(k\right)\left(\widehat{e}\left(k\right)+\frac{\Delta \widehat{e}\left(k\right)}{2}\right)\hfill \\ & ={\left(\frac{\partial \widehat{e}\left(k\right)}{\partial W_1}\right)}^T\Delta W_1\left[\widehat{e}\left(k\right)+\frac{1}{2}{\left(\frac{\partial \widehat{e}\left(k\right)}{\partial W_1}\right)}^T\Delta W_1\right]\hfill \\ & =-{\left(\frac{\partial \widehat{y}\left(k\right)}{\partial W_1}\right)}^T\Delta W_1\left[\widehat{e}\left(k\right)-\frac{1}{2}{\left(\frac{\partial \widehat{y}\left(k\right)}{\partial W_1}\right)}^T\Delta W_1\right]\hfill \\ & =-{\eta }_{w1}{\widehat{e}}^2\left(k\right){∥\frac{\partial \widehat{y}\left(k\right)}{\partial W_1}∥}^2+\frac{1}{2}{\eta }_{w1}^2{\widehat{e}}^2\left(k\right){∥\frac{\partial \widehat{y}\left(k\right)}{\partial W_1}∥}^4\hfill \\ & \le -\frac{1}{2}{\eta }_{w1}\left(2-{\eta }_{w1}{\overline{\mathrm{\Psi }}}_1^2\right){\widehat{e}}^2\left(k\right){∥\frac{\partial \widehat{y}\left(k\right)}{\partial W_1}∥}^2\hfill \end{array}$$

Notice that $\Delta V_1\left(k\right)$ can be negative definite when ${\eta }_{w1}$ satisfies Equation (28). Since $V_1\left(k\right)$ is reducing and non-negative, this function is converging to the constant $V_1^{\infty }\ge 0$, as $k\to \infty$; therefore, $\Delta V_1\left(k\right)\to 0$. It can be inferred that ${\lim }_{k\to \infty }\widehat{e}\left(k\right)=0$. ☐

Based on Equations (18) and (22), the Jacobian information algorithm can be derived as follows:

$$\frac{\partial y\left(k\right)}{\partial u\left(k\right)}\cong \frac{\partial \widehat{y}\left(k\right)}{\partial u\left(k\right)}=\frac{\partial \widehat{y}\left(k\right)}{\partial \widehat{x}\left(k\right)}·\frac{\partial \widehat{x}\left(k\right)}{\partial u\left(k\right)}={\widehat{b}}_0\sum _{m=1}^{p}m{\widehat{\sigma }}_m{\widehat{x}}^{m-1}\left(k\right)$$

(29)

In this part, the presented modified MFACapproach can be used in the SOFC to obtain safe fuel utilization and meet operating constraint conditions under the condition that its current I and voltage $V_{dc}$ can be available.

3.3. Adaptive PID Controller Design with Control Input Constraints

Define tracking error $e\left(k\right)=r\left(k\right)-y\left(k\right)-\zeta$, where $\zeta$ denotes a compensation signal, which can adjust the specified reference $r\left(k\right)$ to guarantee $u\left(k\right)$ in the pre-determined constraints; this is designed later.

According to the basic PID Controller method, the modified PID control algorithm is given, which is suitable for this paper in the following form:

$$\begin{array}{cc}\hfill u_c\left(k\right)=& u(k-1)+k_1{x}_{c1}\left(k\right)+k_2{x}_{c2}\left(k\right)+k_3{x}_{c3}\left(k\right)\hfill \end{array}$$

(30)

where $k_i,i=1,2,3$ denotes the PID parameters, which can be formed as vector $K={[k_1,k_2,k_3]}^T$. Additionally:

$$\begin{array}{c}{x}_{c1}\left(k\right)=e\left(k\right)-e(k-1),x_{c2}\left(k\right)=e\left(k\right)\hfill \\ x_{c3}\left(k\right)=e\left(k\right)-2e(k-1)+e(k-2)\hfill \end{array}$$

Since this PID controller parameters $k_i,i=1,2,3$ are considered as positive or negative, exponential functions are introduced to map the $k_i$, which are defined as:

$$\begin{array}{c}\hfill k_i=\left\{\begin{array}{cc}{\alpha }^{{\omega }_i}\hfill & \mathrm{if}{k}_i\ge 0\hfill \\ -{\alpha }^{{\omega }_i}\hfill & \mathrm{if}{k}_i<0\hfill \end{array}\right.i=1,2,3\end{array}$$

(31)

where $\alpha >1$, ${\omega }_i,i=1,2,3$ are estimated parameters and $\alpha$ is used to change the speed for controller parameters $k_i$.

Due to the input constraints (14), the presented adaptive constrained controller can be written in the following form:

$$u\left(k\right)=\mathrm{Sat}\left\{\left(u(k-1)+\mathrm{Sat}\left\{\left(u_c\left(k\right)-u(k-1)\right),T_s{\dot{q}}_{f\min },T_s{\dot{q}}_{f\max }\right\}\right)\right.,\left.{\overline{q}}_{f\min },{\overline{q}}_{f\max }\right\}$$

(32)

Therein, the $\mathrm{Sat}(·)$ function can be defined as:

$$\begin{array}{c}\hfill \mathrm{Sat}(z,z_{\min },z_{\max })=\left\{\begin{array}{cc}{z}_{\max }\hfill & z\ge z_{\max }\hfill \\ z\hfill & z_{\min }<z<z_{\max }\hfill \\ z_{\min }\hfill & z\le z_{\min }\hfill \end{array}\right.\end{array}$$

where $z_{\max }$ and $z_{\min }$ represent the upper and lower $\mathrm{Sat}(·)$ function, respectively. Since dynamic constraints exist in SOFC, an anti-windup compensator needs to be introduced to adjust the specified trajectory $r\left(k\right)$. This compensation signal $\zeta \left(k\right)$ can be defined in the next form [34]:

$$\begin{array}{c}\hfill \zeta \left(k\right)=\varphi \zeta (k-1)+\frac{\partial y\left(k\right)}{\partial u\left(k\right)}(u_c\left(k\right)-u\left(k\right))\end{array}$$

(33)

Therein, $\varphi <1$.

Remark 2.

Since $0\le \varphi <1$ and we assume $\frac{\partial y\left(k\right)}{\partial u\left(k\right)}$, $u_c\left(k\right)$ are bounded, according to the stability criteria for linear discrete-time systems in [35], the compensation signal $\zeta \left(k\right)$ is UUBfor all k. As we all know, the actuator should drive the controlled system with saturation, when the actuator is unable to totally realize the driving instruction $u_c\left(k\right)$. The proposed anti-windup compensator is to guarantee the actuator operates at saturation (magnitude and rate saturation constraint) by adjusting the reference trajectory when the actuator is unable to realize the driving instruction $u_c\left(k\right)$ completely. The proposed method can ensure the convergence of the parameter update laws.

In the process of the adaptive PID constrained control design, this controller parameter can be adjusted along the negative gradient, and the performance index can be written as:

$$\begin{array}{c}\hfill J_{W_2}=\frac{1}{2}{e}^2\left(k\right)\end{array}$$

(34)

where $W_2={[{\omega }_1,{\omega }_2,{\omega }_3]}^T$. The task obviously becomes an optimization problem, which makes the performance index (34) about weighting parameter ${\omega }_i$, that is:

$$\begin{array}{c}\hfill \frac{\partial J_{W_2}}{\partial {\omega }_i}=-k_i{x}_{ci}\ln \alpha \frac{\partial y\left(k\right)}{\partial u\left(k\right)}e\left(k\right)\end{array}$$

(35)

where $x_{ci},i=1,2,3$ are the defined inputs of the PID controller. Using a back-propagation training algorithm, due to Equation (35), the adaptation law of ${\omega }_i$ is derived as:

$$\begin{array}{cc}\hfill {\omega }_i(k+1)& ={\omega }_i\left(k\right)+\Delta {\omega }_i\left(k\right)\hfill \\ & ={\omega }_i\left(k\right)-{\eta }_{w2}\frac{\partial J_{W_2}}{\partial {\omega }_i}\hfill \\ & ={\omega }_i\left(k\right)+{\eta }_{w2}{k}_i{x}_{ci}\ln \alpha \frac{\partial y\left(k\right)}{\partial u\left(k\right)}e\left(k\right)\hfill \end{array}$$

(36)

From the update law of (36), we can obtain:

$$\begin{array}{c}\hfill \Delta W_2={\eta }_{w2}e\left(k\right)\frac{\partial y\left(k\right)}{\partial u\left(k\right)}\frac{\partial u_c\left(k\right)}{\partial W_2}\end{array}$$

(37)

where:

$$\begin{array}{c}\hfill \frac{\partial u_c\left(k\right)}{\partial W_2}=\ln \alpha {[k_1{x}_{c1},k_2{x}_{c2},k_3{x}_{c3}]}^T\end{array}$$

Theorem 2.

Define ${\overline{\mathrm{\Psi }}}_2={\max }_k∥\frac{\partial y\left(k\right)}{\partial u\left(k\right)}∥$, ${\overline{\mathrm{\Psi }}}_3={\max }_k∥\frac{\partial u_c\left(k\right)}{\partial W_2}∥$. The tracking error $e\left(k\right)$ converges asymptotically to zero if ${\eta }_{w2}$ is chosen as:

$$\begin{array}{c}\hfill 0<{\eta }_{w2}<\frac{2}{{\overline{\mathrm{\Psi }}}_2^2{\overline{\mathrm{\Psi }}}_3^2}\end{array}$$

(38)

Proof of Theorem 2.

Define a Lyapunov function as $V_2\left(k\right)=\frac{1}{2}{e}^2\left(k\right)$; from (37), we get:

$$\begin{array}{cc}\hfill \Delta V_2\left(k\right)& =V_2(k+1)-V_2\left(k\right)\hfill \\ & =\Delta e\left(k\right)\left(e\left(k\right)+\frac{\Delta e\left(k\right)}{2}\right)\hfill \\ & =-\frac{\partial y}{\partial u}{\left(\frac{\partial u_c\left(k\right)}{\partial W_2}\right)}^T\Delta W_2\times \left[e\left(k\right)-\frac{1}{2}\frac{\partial y}{\partial u}{\left(\frac{\partial u_c\left(k\right)}{\partial W_2}\right)}^T\Delta W_2\right]\hfill \\ & =-{\eta }_{w2}e\left(k\right){∥\frac{\partial y\left(k\right)}{\partial u\left(k\right)}∥}^2{∥\frac{\partial u_c\left(k\right)}{\partial W_2}∥}^2\times \left[e\left(k\right)-\frac{1}{2}{\eta }_{w2}e\left(k\right){∥\frac{\partial y\left(k\right)}{\partial u\left(k\right)}∥}^2{∥\frac{\partial u_c\left(k\right)}{\partial W_2}∥}^2\right]\hfill \\ & =-{\eta }_{w2}{e}^2\left(k\right){∥\frac{\partial y\left(k\right)}{\partial u\left(k\right)}∥}^2{∥\frac{\partial u_c\left(k\right)}{\partial W_2}∥}^2+\frac{1}{2}{\eta }_{w2}^2{e}^2\left(k\right){∥\frac{\partial y\left(k\right)}{\partial u\left(k\right)}∥}^4{∥\frac{\partial u_c\left(k\right)}{\partial W_2}∥}^4\hfill \\ & \le -{\eta }_{w2}\left(1-\frac{{\eta }_{w2}{\overline{\mathrm{\Psi }}}_2^2{\overline{\mathrm{\Psi }}}_3^2}{2}\right)\times e^2\left(k\right){∥\frac{\partial y\left(k\right)}{\partial u\left(k\right)}∥}^2{∥\frac{\partial u_c\left(k\right)}{\partial W_2}∥}^2\hfill \end{array}$$

We can find that $\Delta V_2\left(k\right)$ can be negative definite in the variables $e\left(k\right)$ when ${\eta }_{w2}$ satisfies Equation (38). Since $V_2\left(k\right)$ is a decreasing and non-negative function, the function is converging to a constant of zero, i.e., $V_2^{\infty }\ge 0$, as $k\to \infty$; therefore, $\Delta V_2\left(k\right)\to 0$. We can infer that ${\lim }_{k\to \infty }e\left(k\right)=0$.

To provide a direct view impression of the presented control strategy, the structure chart is shown in Figure 5.

4. Simulation Results

In this section, the proposed adaptive constrained control method based on the Wiener-type neural network presented in the above sections is used to obtain safe fuel utilization and make the SOFC system satisfy the operating constraints under the condition of changing load I. At the same time, the voltage output $V_{dc}$ is guaranteed to be steady.

The effective and safe power source will provide stable electric energy for plants; thus, we must guarantee that output voltage $V_{dc}$ is the desired constant value as much as possible, while the external current load may have a negative impact on $V_{dc}$ for the SOFC for some time. In conventional working conditions, the current load I is set as 300 A, and our expected stable electric tension $V_{dc}$ can be set as 332.8 V. The desired fuel utilization is $\left[{\rho }_{\min },{\rho }_{\max }\right]=[0.7,0.9]$ [20]. The maintained operating constraints are magnitude saturation $\left[q_{f\min },q_{f\max }\right]=[0.6,1.2]\mathrm{mol}/\mathrm{s}$ and rate saturation $\left[{\dot{q}}_{f\min },{\dot{q}}_{f\max }\right]=[-0.7,0.7]\mathrm{mol}/{\mathrm{s}}^2$ [20].

For the simulations, the sample time is $T=1s$, and we choose orders of $n_a=2$, $n_b=1$, $n_c=2$ and $p=2$ for the Wiener-type neural network. Notice that open-loop input-output data sequences can obtain the initial values $a_i\left(0\right)$, $b_j\left(0\right)$, $c_l\left(0\right)$ and ${\sigma }_m\left(0\right)$. Moreover, open-loop input-output tested data samples are achieved by exciting the open-loop SOFC control system using given sinusoidal signals $0.7823+0.3\sin \left(0.5t\right)\sin t$ and $300+50\sin \left(0.03t\right)\sin \left(0.04t\right)$ for the fuel and the current demand (i.e., $q_f$ and I), respectively. $a_i\left(0\right)$, $b_j\left(0\right)$, $c_l\left(0\right)$ and ${\sigma }_m\left(0\right)$ can be obtained as:

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{a}_1& a_2\\ b_0& b_1\\ c_1& c_2\\ {\sigma }_1& {\sigma }_2\end{array}\right]=\left[\begin{array}{cc}-0.4573& -0.2286\\ -0.0018& 0.0018\\ 0.1278& 0.1278\\ 1.6567& -0.6735\end{array}\right]\end{array}$$

For the proposed control law, $\alpha$ is chosen as mathematical constant $\mathrm{e}$. The learning rates are selected as ${\eta }_{w1}={[0.7,0.7,0.8,0.8]}^T$ and ${\eta }_{w2}={[15,10,0.1]}^T$. The initial values of the proposed controller are the same as the traditional PID parameters $K_p,K_i,K_d$, i.e., $k_1\left(0\right)=K_p=10$, $k_2\left(0\right)=K_i=2$ and $k_3\left(0\right)=K_d=0.01$.

The simulation results are given from Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10. Firstly, the load current is assumed to change as the curve in Figure 6A with three step signals combined. Under the most stringent step signal conditions, the proposed control algorithm presents a favorable control performance, which is shown in Figure 6B. It shows that once there exists a large mutation of the load current, the corresponding oscillations occur on the stack output voltage, yet disappear soon under the proposed control method, and the voltage tends to be stable. The curves of the natural gas flow rate (i.e., the input) and the acceleration are drawn in Figure 7. As can be seen in Figure 7A, the input magnitude constraint generates a saturation phenomenon around Step 200, and the acceleration is within the confining scope. Once the constraints occur, the proposed anti-windup compensator takes effect, as shown in Figure 8A, at the same time to compensate the performance. The Jacobian information and fuel utilization are depicted in Figure 8B,C, which indicates that the fuel utilization is satisfactory within the excellent range under the proposed control strategy. The curves of the exponential PID parameters and the Wiener-type neural network weight coefficients are presented in Figure 9 and Figure 10, respectively.

Moreover, the conventional PID controller is adopted in the simulation in contrast with the proposed control method. The comparisons focus on the following points:

• In Figure 6B, it can be vividly seen that the voltage under PID control has a greater oscillation amplitude and slower convergence speed than the proposed control method.
• In Figure 7A, the natural gas flow rate under PID control is not as stable as the one under the proposed control method when facing the sudden load disturbance.
• In Figure 8C, the fuel utilization under PID control is beyond the permitted scope, which will have a bad effect on the cell life, while the one under the proposed control method is within an appropriate scope.

To summarize, the proposed adaptive tracking constrained control strategy based on the Wiener-type neural network possesses better control performance on the investigated SOFC system in contrast with the PID control with a smaller oscillation amplitude and faster convergence speed and can be applied in the control of nonlinear systems with a complicated dynamics well.

5. Conclusions

In this paper, a new adaptive tracking constrained control strategy based on the online-updated Wiener-type neural network has been proposed for the SOFC, in which we provide a new effective alternative to overcome the control difficulty of the SOFC. Firstly, the nonlinear dynamical model of SOFC is investigated, and the volt-ampere characteristics for the SOFC are presented. Secondly, the dynamic model of the SOFC is approached using the Wiener model formulation, and the parameters of the model are estimated via adaptive Wiener model identification. Thirdly, an adaptive exponential PID controller is designed to achieve the tracking performance of voltage, with an anti-windup compensator utilized to eliminate the harm of saturation considering the input magnitude and rate constraints. Fourthly, the stability of the proposed control system is proven via the Lyapunov function. Finally, the simulation is carried out in the MATLAB/Simulink environment, and the simulation results are analyzed to show that the proposed adaptive tracking constrained control approach possesses good effectiveness when dealing with nonlinear systems with complicated dynamics. Future work will include solving the problem of output constraints, improving our control algorithm to obtain a better tracking performance, and so on.