Improved Model Predictive Control for the Flow Field in a Wind Tunnel with Model Uncertainty

Improved Model Predictive Control for the Flow Field in a Wind Tunnel with Model Uncertainty

Journal Pre-proof

Improved Model Predictive Control for the Flow Field in a Wind Tunnel with Model Uncertainty Yanming Liu, Junjie Wang, Xiaotao Liu, Xiaoping Li, Bo Yao, Lihao Song PII: DOI: Reference:

S0016-0032(19)30891-9 https://doi.org/10.1016/j.jfranklin.2019.12.004 FI 4321

To appear in:

Journal of the Franklin Institute

Received date: Revised date: Accepted date:

10 June 2019 13 October 2019 1 December 2019

Please cite this article as: Yanming Liu, Junjie Wang, Xiaotao Liu, Xiaoping Li, Bo Yao, Lihao Song, Improved Model Predictive Control for the Flow Field in a Wind Tunnel with Model Uncertainty, Journal of the Franklin Institute (2019), doi: https://doi.org/10.1016/j.jfranklin.2019.12.004

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd on behalf of The Franklin Institute.

Improved Model Predictive Control for the Flow Field in a Wind Tunnel with Model Uncertainty✩ Yanming Liua , Junjie Wanga , Xiaotao Liua,∗, Xiaoping Lia , Bo Yaoa , Lihao Songa a Key Laboratory of Information and Structure Efficiency in Extreme Environment, the Ministry of Education of China, and School of Aerospace Science and Technology, Xidian University, Xian, 710071, China

Abstract This paper designs an improved model predictive control (MPC) strategy for the flow field in a transonic wind tunnel to quickly mitigate the overshoot and smooth the input and output signals under model uncertainty. A flow field control model is established first, followed by model discretization and the general MPC structure. Then, an MPC robust tuning approach is presented, which includes a bending temporal matrix design, and an adaptive robust temporal tuning algorithm. The bending temporal matrix is proposed to reduce the variability in the input profile in the temporal domain by penalizing high-frequency components, thus indirectly smoothing the input and output trajectories. The proposed adaptive robust temporal tuning algorithm can obtain an optimal tuning ratio by maximizing the control performance index of the system under model uncertainty. Simulation results show that the improved MPC controller yields better performance than a traditional MPC controller in mitigating the overshoot and smoothing the input profiles. Keywords: Model predictive control, Wind tunnel, Bending temporal matrix, Model uncertainty, Robust tuning. 2010 MSC: 00-01, 99-00 ∗ Corresponding

author Email addresses: [email protected] (Yanming Liu), [email protected] (Junjie Wang), [email protected] (Xiaotao Liu), [email protected] (Xiaoping Li), [email protected] (Bo Yao), [email protected] (Lihao Song)

Preprint submitted to Journal of the Franklin Institute

January 8, 2020

1. Introduction A transonic wind tunnel is a large tube with air moving inside, which is used for aerodynamic testing scale models of aircraft and spacecraft, in the speed region from zero to Mach 1.3 [1, 2]. Aerodynamic testing data of scale 5

models are measured at a given Mach number, stagnation pressure, temperature and angle of attack. The stagnation pressure (P0 ) and Mach number (M ) are expected to remain stable in a predefined state during transonic wind tunnel tests [3]. Therefore, it is important to keep the Mach number constant at a predefined set point, because most measured variables vary with the Mach

10

number for the given test conditions of stagnation pressure and temperature. Without some kind of control, the variation in the Mach number could be as high as 8% [4]. The proportional-integral-derivative (PID) controller was originally developed as an effective method for flow field control in the transonic wind tunnel

15

process [5]. However, the effectiveness of the PID controller is not very satisfactory when dealing with the perturbation of the flow field caused by a variation in the attack angle [2, 3]. Wind tunnel control systems have had the characteristics of a fast response and good steady performance since Tang [6] optimized the error weighting coefficient matrix by a genetic algorithm. Wang [7] proposed

20

a direct adaptive decoupling controller for the wind tunnel control process to reduce the computational burden and improve the system performance. Fu [8] developed a neural network (NN) adaptive decoupling control approach for the wind tunnel process. Simulation results show that the NN controller yields better performance and effectiveness than the PID controller. However, the dis-

25

advantages of [6, 7, 8] are that quite a large amount of data are required and consumed for the wind tunnel process. On the other hand, MPC has been widely applied to industrial processes, especially for the multivariable control systems [9, 10]. With the available computing power of modern computers, solving complex dynamic optimization prob2

30

lems online is becoming an increasingly viable option to improve the dynamic performance of process operations [11]. Furthermore, MPC can account for the process constraints and multivariable interactions in the optimization problem, such as multi-input multi-output (MIMO) systems. In summary, MPC has the ability of handling multiple variables for the flow field in wind tunnel. Ronald

35

[2] first applied the MPC strategy to a wind tunnel and proved that an MPC controller yields better performance than a PID controller. Simulation and experimental results show that the MPC controller yielded an overall performance improvement of 30-60% compared with that of the PID controller, which is normally used for controlling the Mach number. Zhang [3] applied MPC

40

with a feed-forward strategy for an angle of attack (MPC-FFA) approach to the transonic wind tunnel. Experimental tests demonstrate that the proposed MPC-FFA approach had a higher control precision and a lower test cost than the conventional MPC and PID approaches. However, excessive overshoot leads to jumping parameters and internal flow

45

field disturbance for the wind tunnel, and actuation at high frequencies brings about a decline in the control performance [13, 16], which are not concerned in the above literature. In addition, modeling errors are unavoidable due to the nonlinear variation of the physical characteristics of the wind tunnel, thus the model uncertainty problem occurs. Aiming at the above problems, an improved

50

MPC control strategy for the wind tunnel control process is proposed in this paper. The main contributions of this paper are as follows: (1) A bending temporal matrix for the MPC controller is designed to smooth the input and output profiles, which takes into account the past and predicted control horizon input vectors. The newly designed bending temporal matrix

55

penalizes the high-frequency components more in the frequency domain so that the variability in the input profile is reduced in the temporal domain. (2) A novel adaptive robust temporal tuning algorithm is developed to obtain the optimal robust tuning ratio. The proposed tuning algorithm can obtain the optimal tuning ratio by maximizing the control performance under

3

60

model uncertainty. The remainder of this paper is organized as follows. Section 2 introduces the flow field model and the MPC controller. Then, an MPC robust tuning algorithm is proposed in Section 3. Section 4 shows the simulation results and analyzes the performance of the proposed MPC algorithm. Finally, Section 5

65

concludes the paper.

2. System description and problem formulation In this section, we first introduce the flow field control model and then describe the model uncertainty and the robust stability condition in wind tunnel. Section 2.1 formulates a flow field control model in wind tunnel. In Section 2.2, 70

an MPC controller for the wind tunnel is constructed. Section 2.3 formulates the model uncertainty problem. Finally, the robust stability analysis of the flow field control system is studied in Section 2.4. 2.1. Flow field control model in the wind tunnel The research objectives are Mach number control and stagnation pressure control in the transonic wind tunnel test section [3]. The stagnation pressure regulator P0 remains 1.5 × 105 Pa , and the Mach number value changes from 0 M to 1.2 M with a 0.1 M incremental. Considering the high dynamic and nonlinear characteristics of the transonic wind tunnel, the strategy of multiple models can be applied to the wind tunnel to cover the range of all possible variations of the control models. A simplified block diagram of the flow field control system for the wind tunnel is shown in Fig. 1. A continuous model in the s-domain for the flow field in the transonic wind tunnel could be established as follows      y1 (s) G11 (s) G12 (s) u1 (s)  =  . y2 (s) G21 (s) G22 (s) u2 (s)

(1)

Without loss of generality, the angle of attack remains zero degree, and

the Mach number ranges from 1.1 M to 1.2 M are chosen as a representative 4

Fig. 1: The block diagram of the flow field control system in wind tunnel.

research objective. Finally, the nominal model is given by      20 2.6 −0.4s −0.4s e e u (s) y1 (s) 2 1 (2.5s+1) .  =  5s+1   −13s+6.5 −0.4s 5.8 −0.4s e e u (s) y2 (s) 2 2 (2.8s+1) 3.3s+1

(2)

It is necessary to discretize the continuous model of (1) to apply the MPC

strategy in the wind tunnel [14]. The zero-order hold method is adopted to discretize (1) with sampling periods of Ts = 0.13 s. Then, the discretized model can be represented as follows: x(k + 1) = Ax(k) + Bu(k),

(3)

y(k) = Cx(k). 2.2. MPC Controller 75

For further analysis, an MPC controller for the flow field control system is illustrated in this section based on the industrial transonic wind tunnel process. Fig. 2 shows a block diagram of the designed MPC control system. Yref is the

Fig. 2: The block diagram of the MPC control system.

output target vector, which contains the desired Mach number set points and 5

stagnation pressure set points. The activation function Fr is used to smooth 80

the trajectory of Yref , which is illustrated in Section 3.2. Yref represents the output reference vectors of the Mach number and stagnation pressure. The optimization objectives of flow field control in the transonic wind tunnel system are minimizing the deviations from the reference and control effort [3, 15, 16]. Considering the physics of the wind tunnel process and safety, the state and input constraints must be incorporated [17]. To conclude, the MPC optimization problem is posed as follows:

2

minJ(k) = Yˆ − Yref + kU k2Q2 ∆u

Q1

Hp X = [ˆ y (k + j) − yref (k + j)]T Q1 [ˆ y (k + j) − yref (k + j)] j=1

+

H c −1 X

[uT (k + j)Q2 u(k + j)],

j=0

s.t.

(4)

min{Hc ,i}

x ˆ(k + i) = Ai x ˆ(k) +

X j=1

Ai−j B∆u(k + j − 1),

yˆ(k + i) = C x ˆ(k + i), umin ≤ u(k) ≤ umax , ∆umin ≤ ∆u(k) ≤ ∆umax ,

where U = [u(k), u(k + 1), . . . , u(k + Hc − 1)]T , Yˆ = [ˆ y (k + 1), yˆ(k + 2), . . . , yˆ(k + Hp )]T .

(5)

Yˆ , Yref , U , Hp and Hc are the output prediction vector, reference signal vector, input vector, prediction horizon and control horizon, respectively. umin and umax are lower and upper input limits, respectively. ∆umin and ∆umax are the 85

minimum and maximum change between successive control actions, respectively. u(k + j) is the input variable at time k + j. yˆ(k + i) and x ˆ(k + i) are the output predictions and state at time k + i, respectively. Q1 and Q2 are the block diagonal weighting matrices of the appropriate dimensions. Typically, the weighting matrix Q1 affects the control properties, and Q2 affects the behavior

90

and steady-state performance [13]. When U is determined, the first value of U 6

is sent to the MPC controller as the input. Then, the controller repeats the above steps at each circle but shifts by one step ahead. Based on the receding horizon strategy, the optimization problem shifts forward with the quadratic programming problem that is solved at each iteration step. 95

2.3. Model uncertainty analysis Due to the nonlinear variation of external conditions and disturbances, model uncertainty problem is unavoidable for each of segment of the flow filed control models in wind tunnel. The flow field control mathematical model cannot be obtained accurately and is usually identified by the input-output data in a prac-

100

tical scenario. Therefore, the identified mathematical model in (1) is only an approximation of the real flow field control process. Thus, it is inevitable that the identified nominal model Gnom (s) is different from the real process Gp (s) in the wind tunnel process modeling. Model uncertainty problem results in difficulties in predicting the future states of the control process. If not properly

105

handled, these perturbations may destabilize the MPC controller. In view of the model mismatch between Gnom (s) and Gp (s), model uncertainty problem is considered, and the real process from u(k) to y(k) in the s-domain is assumed to have the following form: Gp (s) =

gp (s + c) , (Tp1 s + 1)(Tp2 s + 1)

(6)

where gp represents the gain of the real process, and Tp1 and Tp2 are the time constants of the real process. Typically, Tp2 = 0 for the first-order model. Then the possible perturbed parameters of the plant model could be denoted such that gp ∈ [(1 + η)gnom , (1 + η)gnom ], (7)

Tp1 ∈ [(1 + η)Tnom1 , (1 + η)Tnom1 ], Tp2 ∈ [(1 + η)Tnom2 , (1 + η)Tnom2 ],

where η, η, and gnom are the minimum uncertainty level value, maximum uncertainty level value, and nominal process gain, respectively. Tnom1 and Tnom2 7

are the time constants of the nominal model. The model uncertainty region ζ is defined as [η, η], which means that the wind tunnel control plant model parameters can be chosen from the above intervals. Given the perturbed parameters of the plant model defined in (7), the model uncertainty can be considered as the set of possible perturbations in the nominal flow field control model. Then, the real flow field control process in the z-domain can be denoted as: (8)

Gp (z) = Gnom (z) + ∆G(z), and 

∆G(z) = 

∆G11 (z)

∆G12 (z)

∆G21 (z)

∆G22 (z)



,

(9)

where Gp (z) represents the real process, and Gnom (z) is the nominal system. Note (8), for all Gp (z) under model uncertainty, i.e., Gp (z) ∈ Π. ∆G(z) is constructed as the additive uncertainty of the subsystems, which is shown in (9). The real process model uncertainty index is defined as δ, which is the max-

110

imum singular eigenvalue of the subprocess model uncertainty in the frequency domain. Typically, for the flow field control process model, the index of the model uncertainty with model uncertainty level ζ and dynamic frequencies ω is shown in Fig. 3.

Fig. 3: δ with ζ and ω.

8

2.4. Robust stability analysis 115

A successfully designed control system should be always able to maintain stability and good performance in the presence of model uncertainty. The uncertainty problem requires the MPC controller to have better robustness to adapt to the model mismatch. The robust stability (RS) condition ensures that the performance of the closed-loop system lies within the specified bounds if the

120

real process is different from its nominal model, even in the worst case condition. Considering the model uncertainty, the closed-loop wind tunnel control system is robustly stable for all Gp (z) in the z-domain if the following RS condition is satisfied:



Tud (z)∆G(z)

∞

(10)

< 1,

where Tud (z) is the control sensitivity function [13], and ∆G(z) is the additive uncertainty depicted in (8). In particular, for the flow field control system in the wind tunnel, Tud (z) can be depicted as follows   Tud11 (z) Tud12 (z) , Tud (z) =  Tud21 (z) Tud22 (z)

(11)

where Tud11 (z), Tud12 (z), Tud21 (z), and Tud22 (z) are the sensitivity functions

Fig. 4: δ with γ and ω.

of the transfer functions [13]. The closed-loop system sensitive function value is defined as α, which is the maximum eigenvalue of a subsystem’s sensitive function in the frequency domain. Let Q2 = q2 I and Q1 = q1 I, where q2 and 9

q1 are the scalar weight values, and I is the adaptive identity matrix. γ is the tuning ratio, where γ = q2 /q1 . The closed-loop system sensitive function value α with γ and ω is shown in Fig.4. The intuitive understanding for the MPC controller is that γ affects the dynamical behavior and steady-state performance. Then, the RS condition in (10) in the frequency domain for the flow field control system can be reformulated as max sup(|Tud (γ, ej2πω )|) < (max |∆G(ζ, ej2πω )|)−1 , ∀ω. γ

ω

ω

(12)

Generally, the smaller the sensitive function value of the closed-loop system is, the better the robustness performance, and vice versa. However, according to the RS condition in (12), the maximum system sensitive function value α should be strictly smaller than δ.

125

3. Robust tuning algorithm In this section, a novel MPC robust tuning algorithm is presented. In Section 3.1, a bending temporal matrix Bt is designed to reduce the high frequency component in the temporal frequency domain. In Section 3.2, an activation function Fr is designed to optimize the output target step signal. Section 3.3

130

introduces the adaptive robust temporal tuning algorithm. 3.1. Automatic bending temporal matrix design To smooth the input trajectory, two modifications are made to the cost function associated with Q2 . First, the input profile term in the past control horizon is added to the cost function as a reference item. Second, a temporal

135

frequency tuning method is adopted for the cost function associated with Q2 . In (4), the cost function associated with Q2 is u(k + j)T Q2 u(k + j) = [Sb u(k + j)]T q2 [Sb u(k + j)],

(13)

where Q2 = SbT q2 Sb , and Sb is the “bending moment matrix” [18, 19]. For the input profile, the first and second order derivatives are incorporated into

10

Sb [20], and the bending behavior is penalized in the cost function. He [16] penalized the high frequency components by constructing Sb in the spatial do140

main. As Sb is a rectangular circulant matrix, it is expected to penalize the high frequency components by constructing a new Sb in the temporal domain. On the basis of spatial frequency tuning theory, a temporal-frequency tuning method is proposed in this paper. The new Sb matrix is defined as Bt , namely, a bending temporal matrix in the temporal-frequency domain. To conclude, the

145

cost function with the input profile term can be realized by Algorithm 1. Algorithm 1 The improved input term cost function with the Bt matrix. 1: procedure input: Hc , q2 , Sb , Bt 2:

Update the input vectors Upast , Upre

3:

if Nsteps ≤ Hc then

do kSb U kQ2 = (Sb Upre )T q2 (Sb Upre )

4: 5:

end if

6:

if Nsteps > Hc then do kBt Utotal kQ2 = (Bt [Upast , Upre ])T q2 (Bt [Upast , Upre ])

7:

end if

8: 9:

end procedure

In Algorithm 1, Upast = Upre Utotal

h

u(k − Hc ), u(k − Hc + 1), . . . , u(k − 1) h iT = u(k), u(k + 1), . . . , u(k + Hc − 1) , h i = Upast , Upre ,

iT

, (14)

where Upast , Upre , and Utotal , and Nsteps are the input vector in the past control horizon, input vector in the predicted control horizon, input vector in the past and predicted control horizon, and iteration times, respectively. Now, the problem turns to the design of Bt for Algorithm 1. 150

It is necessary to analyze the behavior of the existing Sb , because the method of designing Bt is based on Sb . He [16] obtained the spatial frequency response

11

of the given processes and performed an inverse Fourier transform to obtain the bending spatial circulant matrix Sb . Therefore, when Q2 penalizes the input picketing, more weights are placed on the high frequency components of the 155

input profile, which works as a low pass filter to suppress the high frequency components in the input profile [12]. Thus, He penalized the high frequency components by constructing the bending spatial matrix Sb , whose stop band can be adjusted according to different spatial responses in the spatial domain. Based on the method of designing the bending spatial matrix Sb in the spatial

160

frequency domain, a new approach is proposed to design Bt in the temporal frequency domain, whose stop band can be adjusted according to different subsystems to achieve better performance [12, 13, 22]. The method to generate the bending temporal matrix Bt is summarized in Algorithm 2. Fig. 5(a) shows the Algorithm 2 The method of designing the Bt matrix 1: procedure input: G11 , G12 , G21 , G22 2:

Calculate the normalized spectrum response by numerical methods [16];

3:

Calculate low temporal frequency components based on the lowest normalized frequency response;

4:

Design the high pass filter as the response of Bt with preset a stop band.

5:

Process the obtained spectrum content.

6:

Conduct an inverse Fourier transform of the spectrum sequence obtained in line 5 as the weighting matrix Bt .

7:

end procedure

normalized frequency responses of the subsystems and the frequency responses 165

of Bt . Bt can be obtained based on Algorithm 2, and the structure of Bt is shown in Fig. 5(b). To conclude, Bt penalizes the high frequency components more so that the variability in the input profiles is suppressed. As Bt is designed appropriately based on the frequency properties of the subsystems, it is expected to provide better performance than Sb . To test the

170

effectiveness of Sb and Bt in terms of system performance, an sensitivity function Tyd (z) is introduced [13]. Tyd (z) is used to analyze the behavior of the 12

4 2 0 −2 10 (a) The normalized subsystem responses

5

0 0

5

10

(b) The structure of Bt

Fig. 5: The subsystem responses and the structure of Bt .

closed-loop system. In Fig. 6, the performances of Sb and Bt are compared though the values of Tyd (z). The sensitivity functions of G11 , G12 , G21 and G22 are presented in Fig. 6(a), Fig. 6(b), Fig. 6(c) and Fig. 6(d), respectively. 175

On one hand, the values of the sensitivity function with Bt have larger gains at high frequencies, indicating that the high frequency components of the input profile can be penalized more. On the other hand, compared with Sb , the values of the sensitivity function with Bt have smaller gains at low frequencies, demonstrating that the penalization of the input components at low frequencies

180

can be significantly decreased [16]. Therefore, Bt reduces the variability in the input profile in the temporal domain by increasing the penalty of the high frequency components and reducing the penalty of the low frequency components in the frequency domain. Generally, the variation in the input profile trajectory is mitigated by Bt , thus indirectly smoothing the output trajectory. However,

185

this improvement occurs at the cost of slowing the respond speed of the system. 3.2. Improved step target signal design This section focuses on the design of improved step signal in the temporal domain. In addition to the MPC controller, the activation function Fr is an important part of the closed-loop control system (see Fig. 2), which is used for

190

smoothing the step signal of the output target signal yref (k). yref (k) affects the value of the cost function associated with Q1 and indirectly affects the

13

1.5 2.0

Sb

Sb Bt

|Tyd(ej2

)|

|Tyd(ej2

Bt

1.0

)|

1.5

1.0

0.5

0.5

0.0

0.0 0

1

2

3

4

[

s

5

0

6

1

2

3

4

[

]

(a) Sensitivity function of G11

s

5

6

]

(b) Sensitivity function of G12 2.0

1.5

Sb Sb Bt

1.5

Bt

)|

|Tyd(ej2

|Tyd(ej2

)|

1.0

0.5

1.0

0.5

0.0

0.0 0

1

2

3

[

4 s

5

6

0

]

1

2

3

[

(c) Sensitivity function of G21

4 s

5

6

]

(d) Sensitivity function of G22

Fig. 6: Performance comparison of Sb and Bt based on sensitivity function.

14

value of the input. Generally speaking, if yref (k) is a step signal, the input signal varies rapidly with the cost function associated with Q1 . Considering the trade-off between the different performance requirements, it is expected that the overshoot value is minimum and that the rise time value lies in a certain tolerable region. A hyperbolic tangent activation function Fr is designed to smooth the output target step signal in the temporal domain. Fr with a longer length results in a longer rise time, while Fr with a shorter length still leads to overshoot. For the MPC optimization problem solution at every step, the best circumstance is to take the whole variability of the reference step signal into the optimization problem. Therefore, the length of Fr is set to Hp . The procedure to design Fr is as follows: let r be the arguments of Fr , where r = −ceil(Hp /2) : 1 : ceil(Hp /2). Then Fr could be represented as Fr (r) = tanh(r) =

sinh(r) er − e−r = r . cosh(r) e + e−r

(15) 0

Finally, the smoothed target step signal yref at every step can be calculated as 0

yref (k + j) = [yref (k + 1) − yref (k)]Fr (j − ceil(Hp /2)) + yref (k).

(16)

Fig. 7 illustrates the effectiveness of Fr smoothing the target step signal in the temporal domain. The method of applying Fr in the temporal domain is 1.2

Target Signal Profiles

195

Conventional step signal Improved step signal

1 0.8 0.6 0.4 0.2 0 0

0.2

0.4

Time

0.6

0.8

1

Fig. 7: Improved target signal profile with Fr in the temporal domain.

summarized in Algorithm 3. 15

Algorithm 3 The method of applying Fr in the temporal domain ˆ , Yref , Hp 1: procedure input: Y 2: 3:

[ˆ y (k + j) − yref (k + j)]T Q1 [ˆ y (k + j) − yref (k + j)]

4:

end if

5:

if Yref is step signal then

6: 7: 8:

200

if Yref is not step signal then

0

0

[ˆ y (k + j) − yref (k + j)]T Q1 [ˆ y (k + j) − yref (k + j) + j)] end if end procedure

Based on Algorithm 3, if yref is a step signal, the improved cost function associated with Q1 actually penalizes the difference between the reference signal and the smoothed target signal. Consequently, Fr can effectively mitigate the overshoot and optimize the output trajectory by smoothing the output target step signal. However, this improvement occurs at the cost of increasing the

205

rising time of the MPC controller. 3.3. Adaptive robust temporal tuning algorithm In this part, the robust MPC tuning problem is considered to maximize the control performance index with the RS condition. An uncertainty region identification method is estimated based on the nominal model. Then the opti-

210

mal robust tuning ratio can be obtained according to the RS condition and the control performance. The key idea of the uncertainty level identification method is to utilize the deviation degree between the real system step response and the nominal system step response. In particular, for the flow field control model in the wind tunnel, the uncertainty level equals the maximum value of each uncertainty level of the subsystems. The model uncertainty region is defined as ζ, which means the parameters of the real process fluctuate with the parameters of the nominal model up and down within the uncertainty levels. The maximum uncertainty

16

level η can be calculated as follows:  t=n  P    −  Greal(i,j) (t) Gnom(i,j) (t)    η = max t=1 , t=n P  (i,j)      Gnom(i,j) (t)

(17)

t=1

where Greal (i,j) and Gnom (i,j) represent the step responses of the real subsystems and the nominal subsystems, respectively. To conclude, the strategy of the uncertainty level estimation method is summarized in Algorithm 4. Fig. 8 shows the effectiveness of Algorithm 4. Accordingly, the tuning problem can be

Algorithm 4 Uncertainty region identify algorithm 1: procedure input: Gnom 2:

Calculate the step responses of Greal and Gnom ;

3:

Calculate η based on (17);

4:

Let η = −η and η = η;

5:

The model uncertainty region ζ is estimated as [η, η];

6:

end procedure

0.7 Real uncertainty level

0.6 0.5

Upper envelope of the estimated uncertainty level Lower envelope of the estimated uncertainty level

Uncertainty level

0.4 0.3 0.2 0.1 0.0 -0.1 -0.2 -0.3 -0.4 5

10

15

20

Number of tests

Fig. 8: The identification results. 215

formulated as: we tune χ by optimizing the γ with the RS condition. χ is the control performance that can be calculated based on Lemma 1. In particular, for the flow field closed-loop system, the closed-loop system control performance index is chosen as the minimum of the subsystem control performance index. 220

The physical meaning of this index is the volume of the truly removed distur17

bance by the controller. A larger χ means a better control performance; i.e., more disturbances are removed. Lemma 1. Let ω be the dynamical frequency. When |Tyd (ej2πω )| first crosses over 1 from 0 to the dynamical Nyquist frequency, the closed-loop system control performance index is defined by

χ=

k=n P k=1

(1 − Tyd (ej2πω ) )dω ω

(18)

,

where dω denotes the frequency resolution. Based on the above analysis, given the real wind tunnel process, the tuning objective can be expressed mathematically as the following constrained optimization problem in (19). max γ

s.t.

min

χ(γ),

Gp ∈Π

(19)

max sup(|Tud (γ, ej2πω )|) < (max |∆G(ζ, ej2πω )|)−1 , ∀ω. γ

ω

ω

To conclude, the adaptive robust tuning algorithm is summarized in Algo225

rithm 5.

18

Algorithm 5 Adaptive robust tuning algorithm 1: procedure input: Gnom , Hp , Hc , bp∗ 2:

Estimate the uncertainty region ζ based on Algorithm 4;

3:

γend ← 0, γstart ← 102

4:

while γstart − γend > bp∗ do

5:

γ ← (γend + γstart ) × 0.5

6:

Calculate Tud (γ) if Tud (γ, ej2πω ) < δ then

7: 8: 9: 10: 11: 12:

Record χ(γ, ζ)

end if γend ← γ if Tud (γ, ej2πω ) ≥ δ then γstart ← γ

13:

end if

14:

γART ← max

15:

end

16:

χrecord (γ)

while

end procedure In Algorithm 5, bp∗ , γstart , γend , and γART are the bisection precision, start

region of the tuning ratio, terminal region of the tuning ratio, and optimal tuning ratio, respectively. Taking Bt into account, Q2 is calculated as Q2 = BtT q1 γBt in Algorithm 5. The proposed adaptive robust tuning algorithm can obtain 230

the optimal robust tuning ratio with different specific model uncertainty levels automatically. Taking Fr and Bt into account, the automatic robust temporal tuning (ARTT) algorithm is concluded in Algorithm 6. Algorithm 6 is an offline tuning strategy that is based on the bisection method, i.e., to effectively decrease the computing time and obtain the best control performance. The effectiveness

235

of Algorithm 6 will be illustrated in the next section with a real flow field control industrial model.

19

Algorithm 6 Automatic robust temporal tuning algorithm (ARTT) 1: procedure input: Hp , Hc , Gnom , Yref , ∆umin , ∆umax , bp∗ 2:

Construct Bt based on Algorithm 2;

3:

Optimize the target value based on Algorithm 3;

4:

Identify the uncertainty region ζ based on Algorithm 4;

5:

Obtain γART based on Algorithm 5;

6:

if Nsteps ≤ Hc then

7:

do kSb U kQ2 = (Sb Upre )T q1 γART (Sb Upre )

8:

end if

9:

if Nsteps > Hc then

10: 11: 12:

do kBt Utotal kQ2 = (Bt [Upast , Upre ])T q1 γART (Bt [Upast , Upre ]) end if end procedure

4. Numerical simulation results and discussion To illustrate the effectiveness of the proposed algorithm, extensive evaluations and comparisons are presented in this section by applying the algorithm 240

to the flow field in the transonic wind tunnel. As ARTT algorithm is based on the He’s [12] user-friendly robust temporal tuning (URTT) approach, it is expected that ARTT tuning approach yields better performance than URTT tuning approach. Different from ARTT approach, only RS condition and Sb matrix are taken into account in URTT approach. The input constraints are

245

umin = 0 and umax = 5, the prediction horizon Hp is 5, and the control horizon Hc is 3. Without loss of generality, the unknown model uncertainty level is set as 0.15. The identified uncertainty region ζ result by Algorithm 4 is 0.15. The computational costs of the proposed tuning method for the wind tunnel control process accounts for 36.674 s on a computer with i7 CPU and 8GB RAM.

250

To compare the performance of URTT approach and ARTT approach, the output and input results for the flow field control model under the same model uncertainty region are shown in Fig. 9, Fig. 10, Fig. 11, and Fig. 12. The output

20

envelope results of the flow field for the wind tunnel under uncertainty level ζ are shown in Fig. 9 and Fig. 10, of which Fig. 9(a) and Fig. 9(b) are for the 255

ach number envelopes, and Fig. 10(a) and Fig. 10(b) are for the stagnation envelopes. As shown in Fig. 9(a) and Fig. 9(b), both the MPC with URTT and the MPC with ARTT can track the target values accurately. Table 1 is prepared on the basis of Fig. 9 and Fig. 10 to quantitatively characterize the controller tracking performance in the temporal domain. The statistical data

260

show that the both of the overshoots and settling times of the outputs of the MPC with ARTT are effectively reduced within the acceptable steady error (0.5%). Compared with the results of the MPC with URTT, the Mach number and stagnation pressure output envelopes of the improved MPC strategy with the ARTT algorithm are smoother. In addition, the overshoot is mitigated

265

effectively by the ARTT algorithm. However, as Fr is applied to the ARTT algorithm to smooth the step signal of the Mach number output target value, the rise time of the Mach number output of the ARTT algorithm becomes a bit longer. 1.3

Nominal

Nominal

Envelope

Envelope Mach number

Mach number

1.3

1.2

1.2

1.1

1.1 0

1

2

3

4

0

5

1

2

3

Time (s)

Time (s)

(a)

(b)

(a) MPC with URTT

4

5

(b) MPC with ARTT

Fig. 9: Mach number value.

The main regulating valve voltage input envelopes and the main exhaust 270

valve voltage input envelopes are shown in Fig. 11 and Fig. 12, respectively. As can be seen, both of the two input envelopes of MPC with ARTT approach are smoother than that of MPC with URTT. In addition, Fig. 11 shows that

21

1.7

1.7

Nominal

Nominal

Stagnation pressure

Stagnation pressure

Envelope

1.6

1.5

1.4

1.3

Envelope

1.6

1.5

1.4

1.3

1.2 0

1.2 0

1

2

3

4

1

2

5

3

4

5

Time (s)

Time (s)

(b)

(a)

(a) MPC with URTT

(b) MPC with ARTT

Fig. 10: Stagnation pressure value.

Table 1: Dynamic performance of URTT [12] and ARTT

Tuning approach

Performance Category

Overshoot

Settling time

Steady error

URTT

Mach number

0.66%

1.82 s

1‰

Stagnation pressure

6.8%

2.34 s

0.7‰

Mach number

0.50%

1.56 s

1.3‰

Stagnation pressure

1.2%

1.30 s

0.0005‰

ARTT

the main regulating valve voltage input envelope region of MPC with ARTT is obviously smaller than that of MPC with URTT. Both of the above improvements arise mainly because of the addition of Bt to the ARTT algorithm. This method penalizes more high-frequency components in the frequency domain, thus indirectly smoothing the input profiles in the temporal domain. 0.5

0.5

Nominal

Nominal Envelope

0.4

Main regulating valve

Main regulating valve

275

0.3

0.2

Envelope

0.4

0.3

0.2

0.1

0.1 0

1

2

3

4

0

5

1

2

3

Time (s)

Time (s)

(a)

(b)

(a) MPC with URTT

(b) MPC with ARTT

Fig. 11: Main regulating valve voltage.

22

4

5

1.4

1.4

Nominal

Nominal 1.3

Envelope Main exhaust valve

Main exhaust valve

1.3

1.2

1.1

1.0

0.9

Envelope

1.2

1.1

1.0

0.9

0.8

0.8 0

1

2

3

4

0

5

1

2

3

Time (s)

Time (s)

(a)

(b)

(a) MPC with URTT

4

5

(b) MPC with ARTT

Fig. 12: Main exhaust valve voltage.

Additionally, Fig. 12 shows that the envelope areas of the main exhaust valve voltage of MPC with ARTT remains almost unchanged compared with those 280

of MPC with ARTT. However, the edges of the main exhaust valve envelope become smooth. This is because that the Mach number is mainly affected by the main regulating valve, not the main exhaust valve.

5. Conclusion This paper develops an improved MPC control scheme for the flow field con285

trol in a transonic wind tunnel, including an MPC controller design and an MPC robust tuning algorithm. The optimal tuning ratio with the identified model uncertainty level can be obtained by the proposed ARTT algorithm. Moreover, ARTT algorithm could mitigate the overshoot and optimize the input and output profiles significantly. The results clearly indicate that the improved MPC

290

with the ARTT algorithm outperforms that if MPC with URTT in overshoot mitigation of the output profiles and high-frequency component reduction of the input profiles. However, experiments for the wind tunnel control system should be conducted to validate the performance of the proposed MPC control scheme. Besides, the limitations of the proposed ARTT algorithm are that the

295

identification accuracy can be challenging and that the computation should be reduced. These extensions will be investigated in future work.

23

Acknowledgment This work was supported by the National Natural Science Foundation of China (grant Nos. 61627901, 61703320, and 61431010), and Doctoral Students’ 300

Short Term Study Abroad Scholarship Fund.

References References [1] W. Yu, and G. Zhang, “Modelling and controller design for 2.4 m injector powered transonic wind tunnel,” presented at the American Control 305

Conference. IEEE, USA, June 6-6, 1997. [2] R.A.M. Soeterboek et al., “A predictive controller for the Mach number in a transonic wind tunnel,” IEEE Control Systems Magazine , vol.11, no. 1, pp. 63–72, Jan. 1991. [3] J. Zhang, et al., “Model Predictive Control for the Flow Field in an In-

310

termittent Transonic Wind Tunnel,” IEEE Transactions on Aerospace and Electronic Systems, vol.54, no. 1, pp. 324–338, Feb. 2018. [4] R. Soeterboek. “Predictive control: a unified approach,” Prentice-Hall, Inc., 1992. pp. 293–306. [5] A.F. Pels. “ Closed-loop Mach number control in a transonic wind-tunnel,”

315

Journal A, vol.30, no.4, pp. 25–32, 1989. [6] Y. Tang et al., “ State space description based predictive control for normal temperature continuous transonic wind tunnel,” 2016 Chinese control and decision conference (CCDC). IEEE, Yinchuan, China, May 28-30, 2016. [7] X. Wang et al., “Multi-model direct adaptive decoupling control with ap-

320

plication to the wind tunnel system,” ISA transactions, vol.44, no.1, pp. 131–143, Jan. 2005.

24

[8] Y. Fu and T. Chai, “Intelligent decoupling control of nonlinear multivariable systems and its application to a wind tunnel system,” IEEE Transactions on control systems technology , vol.17, no.6, pp. 1376–1384, Apr. 325

2009. [9] X. Liu, “Robust model predictive control and distributed model predictive Control: Feasibility and Stability,” Ph.D. dissertation, University of Victoria, CAN, 2016, pp. 2–14. [10] T. Chai et al., “Optimal operational control for complex industrial pro-

330

cesses,” Annual Reviews in Control, vol.38, no.1, pp. 81–92, 2014. [11] E. Matthew et al., “A tutorial review of economic model predictive control methods,” Journal of process control , vol.24, no.8, pp. 1156–1178, Aug. 2014. [12] N. He, “User friendly robust MPC tuning of uncertain paper-making pro-

335

cesses,” Ph.D. dissertation, University of Alberta, CAN, 2017, pp.77–82. [13] J. Fan, “Model predictive control for multiple cross-directional processes: analysis, tuning, and implementation,” Ph.D. dissertation, University of British Columbia, USA, 2003, pp. 81–91. [14] T. Chen, “Optimal sampled-data control systems,” Springer Science &

340

Business Media, 2012, pp.31-49. [15] J. Fan, and G. E. Stewart, “Automated tuning of large-scale multivariable model predictive controllers for spatially-distributed processes,” presented at the Proceedings of the 2006 American Control Conference, Minneapolis, Minnesota, USA, June 14-16, 2006.

345

[16] N. He, D. Shi, F. Michael, B. Johan, and T. Chen, “Robust tuning for machine-directional predictive control of MIMO paper-making processes,” Control Engineering Practice, vol.55, no. 8, pp. 1–12, Oct. 2016.

25

[17] X. Liu et al., “Multistage suboptimal model predictive control with improved computational efficiency,” Journal of Dynamic Systems, Measure350

ment, and Control, vol.136, no.3, pp. 031026, Mar. 2014. [18] R.A. Bartlett et al., “Quadratic programming algorithms for large-scale model predictive control,” Journal of Process Control, vol.12, no.7, pp. 775–795, Oct. 2002. [19] J.G. Vanantwerp et al., “Fast model predictive control of sheet and film

355

processes,” IEEE Transactions on Control Systems Technology , vol.8, no.3, pp. 408–417, May 2000. [20] J. Fan et al., “Approximate steady-state performance prediction of largescale constrained model predictive control systems,” vol.13, no.6, pp. 884– 895, Nov. 2005.

360

[21] J. Fan, “ Two-dimensional frequency analysis of structured uncertainty for multiple array paper machine cross-Directional processes,” Proceedings of the 44th IEEE Conference on Decision and Control , Seville, Spain, Dec. 2005. [22] G.E. Stewart, “Feedback controller design for a spatially distributed sys-

365

tem: The paper machine problem,” IEEE Transactions on Control Systems Technology , vol.11, no.5, pp. 612–628, Sept. 2003.

26