Robust Tuning of Machine Directional Predictive Control of Paper

Mar 19, 2015 - On the basis of the visualization technique and the unmodality/monotonicity properties of the performance indices with respect to the t...
0 downloads 4 Views 4MB Size
Article pubs.acs.org/IECR

Robust Tuning of Machine Directional Predictive Control of Paper Machines Dawei Shi,*,†,§ Jiadong Wang,† Michael Forbes,‡ Johan Backström,‡ and Tongwen Chen† †

Department of Electrical and Computer Engineering, University of Alberta, Edmonton, Alberta T6G 2V4, Canada Honeywell Process Solutions, 500 Brooksbank Avenue, North Vancouver, British Columbia V7J 3S4, Canada



ABSTRACT: In this work, a parameter tuning problem of two-degrees-of-freedom model predictive control of industrial papermaking processes is explored to achieve satisfactory time-domain robust closed-loop performance in terms of worst-case overshoots and worst-case settling times, under user-specified parametric uncertainties. An efficient visualization method is first developed to characterize the set of time-domain closed-loop responses in the presence of parametric model−plant mismatch. On the basis of the visualization technique and the unmodality/monotonicity properties of the performance indices with respect to the tuning parameters, the feasibility of the tuning problem can be analyzed, and a three-step iterative line-search based automatic tuning algorithm is proposed to determine the controller parameters that meet the time-domain performance requirements robustly for the given parametric uncertainty specifications. The effectiveness of the algorithm is illustrated by applying the results to a process from stock to conditioned weight in an industrial paper machine and by comparing the performance of the algorithm with that of brutal search.

1. INTRODUCTION Model predictive control (MPC) has been widely applied for quality control of paper-making processes.2,3 Different from the PID control, another well-known control strategy that is normally applied at the base level of a plant, the MPC controllers are typically used in the supervisory mode at the high level,4 for example, in the case of paper-machine quality control systems. Despite its vast applications, a significant gap exists between the theory of MPC and its practical usage,5−7 due to the underlying mathematical difficulties in addressing the properties (e.g., stability, feasibility, robustness) of the approach. Among others, a crucial problem that the industrial MPC practitioners are eager to solve is MPC parameter tuning for desired performance in terms of well-known indexes (e.g., overshoots and settling times), which is the focus of the work in this paper. The tuning problem considered in our work arisen in the machine directional (MD) MPC (which we refer to as MDMPC hereafter) used for quality control of paper machines in the pulp and paper industry; however, similar problems should exist for MPC applications in other industrial processes, for example, the oil and gas industry,8,9 due to the MPC formulation and its extensive applications. The motivation of this work mainly stems from the requirements from the end users of the commercial MPC design and tuning tools in the pulp and paper industry. The practitioners/ operators of paper machines are normally more familiar with the classical performance indices on step responses, for example, the overshoot, rising time, and settling time. Although the expressions of these indices are well known for standard secondorder systems, the expressions of them for the closed-loop system controlled by MPC are difficult to find and are too complicated to be considered for controller tuning. Moreover, the quantitative relationship of these performance indices with the indices that are easier to work with in modern control theory (e.g., H2/H∞ norms) does not seem to exist. Therefore, when the © XXXX American Chemical Society

performance specifications are provided in terms of overshoots and settling times, the tuning problems cannot be solved by extending the results in H2/H∞ theory, although stability margins can still be obtained based on H∞ analysis (as will be shown in this work). Consequently, the structure of the resultant tuning problem becomes implicit and mathematically difficult to handle. The end users also require the robustness of tuned MPC controllers, so that satisfactory responses (in terms of “worstcase” overshoots and “worst-case” settling times) can be maintained when model mismatch or uncertainty occurs during the process operation, which invariably happens due to the modeling errors and the change of operating conditions. In particular, the end users are much more comfortable with the intuitive specifications of uncertainty for the parameters of a first order plus dead time (FOPDT) model than with unstructured uncertainty, which is more extensively studied in the robust control literature;10 for the former type of uncertainty, a lot of problems considered are unfortunately NP-hard.11 This further adds to the difficulty in designing the tuning tools, as the worstcase overshoot and worst-case settling time need to be identified and optimized in the tuning procedure. In general, MPC parameter tuning problems are essentially difficult to solve, even for the unconstrained case. As is shown in Table 1, a number of interesting attempts have been reported in the literature.12−18 An in-depth review of the results on MPC tuning was provided in Garriga and Soroush.19 From Table 1, a variety of tuning objectives have been considered for MPC tuning; however, it is still not clear how to automatically determine the MPC parameters so that the user-friendly specifications Received: December 30, 2014 Revised: March 3, 2015 Accepted: March 19, 2015

A

DOI: 10.1021/ie5050583 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Table 1. Comparison of Different Studies on MPC Tuning authors

processes

Al-Ghazzawi et al. (2001)12 Garriga and Soroush (2008)13 Di Cairano and Bemporad (2010)14 Shah and Engell (2010)15 Shah and Engell (2011)16

MIMO linear systems SISO linear systems MIMO linear systems SISO linear systems MIMO linear systems

Toro et al. (2011),17 Vallerio et al. (2014)18

nonlinear systems

design specifications graphical time-domain specifications closed-loop eigenvalue specifications matching the unconstrained MPC with a prespecified controller desired closed-loop poles and zeros closed-loop performance specifications given in terms of a preferred closed-loop transfer function balancing the conflicting terms contained in the MPC cost function

on overshoots, settling times, and parametric uncertainty can be simultaneously fulfilled. Note that although the requirements on overshoots and settling times can be implicitly taken into account using the graphical time-domain specifications in Al-Ghazzawi et al.,12 the issue is that the hard constraint on overshoot cannot be guaranteed and the settling time cannot be directly minimized. On the other hand, one possible way to reduce the difficulty of MPC tuning is to modify the controller structure.20,21 Interestingly, in Bordons and Camacho,20 a generalized predictive controller was introduced for FOPDT processes, and a predictive control strategy for stabilizable linear systems was proposed by Kong et al.;21 the shared benefit of these controller frameworks is that the controller parameters are easier to be adjusted for desired performance requirements. In particular, the controller in Bordons and Camacho20 can be tuned like a PID controller and a single tuning parameter was utilized in Kong et al.21 to adjust the controller performance. On the basis of the aforementioned requirements in the applications of industrial MPC tuning and illuminated by the general idea of modifying the controller structures,20,21 an alternative tuning approach is proposed in this work. The starting point of our approach follows from the two-degree-of-freedom (2-DOF) MD-MPC framework introduced in Chu, Forbes, and Backström.22 Instead of tuning the weighting matrices in the control optimization problem, two filter parameters (which we refer to as λ parameters hereafter) are used to adjust the closedloop performance. We assume the nominal process model is known and consider the parametric uncertainties on the process parameters (namely, the model−plant mismatch), which is straightforward to specify by the end users. The performance requirements are quantified in terms of worst-case overshoots and settling times. As has been discussed earlier, the motivation and target of this work are to meet the needs of technicians and process engineers, rather than to make theoretic advances. Both the nominal process and the real process are assumed to have the FOPDT form, which has been widely used to model machine directional and cross directional dynamics of industrial paper machines.22−24 Rather than the MPC with constraints, an unconstrained MPC is considered in the tuning problem for the following reasons: 1. Time-consuming parameter tuning processes for both offline and online applications lead to unnecessarily high maintenance costs and economic loss on product quality and productivity of industrial manufacturing processes, as the cost for tuning (e.g., on human resources) increases with the amount of tuning time and the manufacturing processes have to stop or run with poorly tuned controllers during the tuning procedure. The computation time of the autotuning algorithm, however, can be reduced by one or two magnitudes if the MPC constraints are ignored, which is of great benefit and importance to industrial tuner design.

2. Experience on paper machine quality control indicates that the performance of the tuned unconstrained MD-MPC efficiently approximates that of the constrained case with the same tuning parameter settings. The intuition is that the constraints in MD-MPC become inactive as the process variables approach their operating points. Thus, the effect of the constraints is ignored here so that we can solely focus on our end users’ needs on specifying parametric uncertainty, resisting unacceptable overshoots, and minimizing settling times. The formulated MPC tuning problem maintains the friendliness of the proposed results to paper-machine operators/end users but further increases the difficulty of the analysis and parameter tuning. Because the inevitable existence of time delays, analytical expressions for the closed-loop responses and the performance indices become too complicated to be used for parameter tuning and in particular, are hard to find. The structure of the resultant underlying optimization problem thus becomes unclear. Moreover, to obtain satisfactory user experience, the computation time for the tuning procedure is very limited. Considering these difficulties, an efficient heuristic approach is utilized to find a solution to the tuning problem. The main results and contributions of this work are summarized as follows: 1. On the basis of the small gain theorem, a sufficient condition for robust stability under the parametric uncertainties is presented, which determines the feasible region of the λ parameters. 2. An efficient method to visualize the time-domain closedloop performance in the presence of the model−plant mismatch specifications is proposed. The efficiency of the technique in terms of computational time and precision is validated via extensive simulations. The technique can also be applied to the visualization of other system variables, for example, the control signal. 3. A three-step iterative algorithm is proposed to automatically determine the λ parameters that meet the timedomain robust tuning specifications based on the specified parametric model−plant mismatch, utilizing the visualization technique and the empirical monotonic/unimodality properties of the performance indices with respect to the λ parameters. 4. The performance of the algorithm is tested extensively for different uncertainty levels and worst-case overshoot specifications on an identified process model that is used for machine directional MPC of a paper machine at an industrial site. For all the examples considered, the algorithm always provides almost-optimal tuning parameters in the sense of constraint satisfaction for the obtained worst-case overshoots and optimality of the resultant worst-case settling time, compared with those obtained by brutal search. B

DOI: 10.1021/ie5050583 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 1. 2-DOF MD-MPC structure.

Letting Δx(k + 1) := x(k + 1) − x(k) and Δu(k + 1) := u(k + 1) − u(k), we construct an augmented state space model as follows:

The tuning method has been designed for quality control of paper machines, but is equally applicable to general industrial processes where end users have limited knowledge of control theory. The proposed algorithms can be extended to the case of MPC with constraints and nonlinear models by incorporating different and more complex process simulators in the performance visualization algorithm (namely, Algorithm 1), but the sacrifice is that the computation time will be increased considerably. In this sense, the general idea of tuning applies not only to the 2-DOF MD-MPC framework, but also potentially to other MPC structures. The rest of the paper is organized as follows: Section 2 presents the problem setup. In Section 3, stability analysis is provided based on the small gain theorem. A performance visualization technique is presented in Section 4. In Section 5, an iterative parameter tuning algorithm is proposed. Section 6 presents the application of the proposed results to a paper machine process in machine directional control of paper machines, followed by some concluding remarks in Section 7.

xa(k + 1) = Aaxa(k) + Ba Δu(k) y(k) = Caxa(k)

where ⎡ Δx(k)⎤ ⎥ xa(k):=⎢ ⎢⎣ Cx(k) ⎥⎦

g e−tds Tps + 1

y(k) = Cx(k)

(3)

x(̂ k + 1) = A 0x(̂ k) + B0 u(k)

(6)

y (̂ k) = C0x(̂ k)

(7)

Let Δx̂(k + 1) := x̂(k + 1) − x̂(k). Similar to the real process, an augmented state space model can be constructed as xâ (k + 1) = Aâ xâ (k) + Bâ Δu(k)

where g, Tp, and td are the process gain, the time constant, and the time delay, respectively. The intuition of this assumption is that this type of process model can be used to effectively approximate the nonlinear system equations around the operating points and the high-order dynamics.25 As state-space models are more convenient in MPC implementation and analysis, we denote the discrete-time state space realization of G(s) as (2)

⎡B⎤ Ba :=⎢ ⎥ ⎣CB ⎦

where g0, Tp0, and td0 are the nominal process gain, the nominal time constant, and the nominal time delay, respectively. Note that this can be done by industrial process identification software, for example, the Honeywell Profit Design Studio. The discretetime state space realization of model 5 can be written as

(1)

x(k + 1) = Ax(k) + Bu(k)

⎡ A 0⎤ Aa :=⎢ ⎣CA 1 ⎥⎦

and Ca := [0 1]. This model will be used in the performance visualization technique proposed in this work. Note that the parameters of the real process are unknown and are subject to change caused by the change of operating points; taking account of this issue, we assume that these parameters lie in a known region determined by the identified model and the user-specified uncertainty levels, as will be introduced later. 2.2. Nominal Model. As the real process G(s) cannot be exactly known, we use G0(s) to denote the FOPDT model identified from the input/output data of the real process: g0 e−td0s G0(s) = Tp0s + 1 (5)

2. SYSTEM DESCRIPTION AND PROBLEM FORMULATION In this section, the 2-DOF controller structure for MD-MPC proposed in Chu, Forbes, and Backström22 is first revisited (see Figure 1), and then the tuning problem is formulated. The closed-loop system is composed of four parts: the real process G(s), the nominal model G0(s), the MPC and the user-specified filters (Fr and Fd), which are introduced in detail in the following. Since overshoots and settling times are defined for set-point changes, we assume d(k) = 0 in this work. 2.1. Real Process. The MD process from u(k) to y(k) is assumed to have the FOPDT form G (s ) =

(4)

y ̂(k) = Câ xâ (k)

(8)

where ⎡ Δx(̂ k) ⎤ ⎥ xâ (k):=⎢ ⎢⎣C0x(̂ k)⎥⎦

⎡ A0 0⎤ ⎥ Aâ :=⎢ ⎢⎣C0A 0 1 ⎥⎦

⎡ B0 ⎤ ⎥ Bâ :=⎢ ⎢⎣C0B0 ⎥⎦

and Ĉ a := [0 1]. This model will be used to get the state predictions in the MPC formulation. 2.3. MPC Solution. For the MPC problem considered, we denote Hu as the control horizon and Hp the prediction horizon C

DOI: 10.1021/ie5050583 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research (1 ≤ Hu ≤ Hp). Based on eq 8, the prediction model can be rewritten as min{Hu, i}

i xâ (k + i) = Aâ xâ (k) +



Fr(z) =

Fd(z) =

bdz −1 1 − adz

−1

z −Td0

T − λ Ts d p0

where ar = e , br = 1 − ar, ad = e , bd = 1 − ad, and Td0 is the discretized version of td0. The state space realizations of Fr(z) and Fd(z) can be readily constructed following a procedure similar to the construction of the prediction model in eq 9. With these filters, the MPC performance can be adjusted by tuning λ and λd and setting Q1 = I, Q2 = Q3 = 0, which simplifies the tuning problem.22 In this formulation, although the MPC cost function does not have any penalty on the control input, which is normally used in the conventional MPC formulations to prevent aggressive control moves, the compromise between the transient performance and the smoothness of the responses can be alternatively adjusted by the λ parameters. This also partially explains how the MPC tuning problem is simplified by the considered framework: the MPC penalizes the states only, and the λ parameters govern the aggressiveness/smoothness of the control moves according to user-specified requirements on overshoots and settling times. 2.5. Tuning Objectives. As the nominal model is identified from input/output data of a real process (see Figure 1), modeling errors are unavoidable. Therefore, we allow end users to specify their belief on the precision of the obtained model by providing the trust regions around the nominal model parameters such that

(9)

where i = 1, 2, ..., Hp. The following cost function is considered: J(Y ̂ , ΔU ) :=(Y ̂ − Yref )T Q 1(Y ̂ − Yref ) + ΔUTQ 2ΔU + (U − Uref )T Q 3(U − Uref )

where Q1, Q2, and Q3 are weighting matrices, ⎤ ⎡ Δu(k) ⎥ ⎢ ⎢ Δu(k + 1) ⎥ ΔU = ⎢ ⎥ ⋮ ⎥ ⎢ ⎢ Δu(k + H − 1)⎥ ⎦ ⎣ u 0 ⋯ 0⎤ ⎥ 1 ⋱ ⋮⎥ ΔU ⋱ ⋱ 0⎥ ⎥ ⋯ 1 1⎦

and Uref and Yref are similarly defined dimension-compatible vectors but are composed of elements of uref and yref, respectively, that provide the references of the corresponding variables at different time instants. For typical commercial MPC implementations, constraints on ΔU and U are also considered for quality control; however, as has been discussed in Section 1, these constraints increase the difficulty of MPC tuning, whereas the numerical evaluation of the constrained quadratic programming problems adds to the computation burden. A detailed comparison on the computation burden for the constrained and unconstrained MPC tuning is presented in Remark 3. To simplify the tuning analysis and to speed up the tuning algorithm, we consider an unconstrained MPC

g ∈ [(1 − r̲ g)g0 , (1 + rg̅ )g0]

(12)

Tp ∈ [(1 − r̲ Tp)Tp0 , (1 + rT̅ p)Tp0]

(13)

td ∈ [(1 − r̲ td)td0 , (1 + rtd̅ )td0]

(14)

where rg ∈ [0,1], rTp ∈ [0,1], rtd ∈ [0,1], rg ∈ [0,∞), rTp ∈ [0,∞), and rtd ∈ [0,∞) are used to describe the levels of the parametric uncertainties. These trust regions will also be denoted as g ∈ [g, g], Tp ∈ [T p, Tp], td ∈ [t d, td] for notational brevity. Note that compared with unstructured uncertainty, the parametric uncertainty is relatively easier to understand and specify for the end users, especially for those without a background on robust control theory. We now present the performance requirements for parameter tuning, and formulate the tuning problem into a constrained optimization problem. Considering the preference of the end users, we quantify the closed-loop performance with the overshoot and settling time. In industrial applications, the requirement on overshoot is that it should not exceed a certain level, whereas a smaller settling time is always more preferable as it quantifies the duration of the transient behavior. In this way, the parameter tuning problem is formulated into the following optimization problem:

minΔU J(Y ̂ , ΔU ) s. t. dynamic equation (9)

z −Td0

T − λTs p0

y (̂ k + i) = Câ xâ (k + i)

⎡1 ⎡1⎤ ⎢ ⎢ ⎥ 1 1⎥ ⎢ U= u(k − 1) + ⎢ ⎢⋮ ⎢⋮⎥ ⎢ ⎢⎣ ⎥⎦ ⎣1 1

1 − ar z

−1

(11)

i−j Aâ Bâ Δu(k + j − 1)

j=1

⎡ y (̂ k + 1) ⎤ ⎢ ⎥ ⎢ y (̂ k + 2) ⎥ ⎥ Ŷ = ⎢ ⋮ ⎢ ⎥ ⎢ ⎥ ⎣ y (̂ k + Hp)⎦

brz −1

(10)

in our industrial MPC tuner design, for which a closed-form solution can be analytically derived by calculating the matrix derivatives and is denoted as Δu(k) = f MPC(Uref, Yref) for notational simplicity. 2.4. Two-DOF Control. In addition to the nominal model and the MPC, two filters Fr and Fd are introduced to aid the parameter tuning.22 These filters are respectively used for filtering the output target, ytgt(k), and the estimated disturbance, d̂(k) := (k) − ŷ(k). With the filtered signals, we calculate the reference trajectory according to

min ts(λ , λd) λ , λd

s. t. ζw(λ , λd) ≤ ζw*

yref (k) = Frytgt (k) − Fdd(̂ k)

(15)

Note that ts(λ, λd) and ζw(λ, λd) here denote the worst-case settling time and worst-case overshoot for g ∈ [g, g], Tp ∈ [Tp, Tp], td ∈ [td, td], which will be further formally defined and evaluated in Section 4. Also, the above formulation automatically rules out the λ parameters that lead to unstable closed-loop

In this way, we have yref(z) = Fr(z)ytgt(z) − Fd(z)d̂(z); Fr(z) and Fd(z) are the so-called reference tracking filter and disturbance rejection filter D

DOI: 10.1021/ie5050583 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 2. Simplified structure of the control loop.

Figure 3. General structure of the closed-loop system.

From the classical results in robust control,26 we know that stability of the closed-loop system, F(N, Δ), is equivalent to the stability of the system in Figure 3b, where M = N11. Therefore, we here examine the M−Δ system instead of the whole system. For the MPC system considered, we have

systems, as a unstable response would imply an unbounded settling time.

3. ROBUST STABILITY ANALYSIS Starting from this section, we discuss the tuning of λ and λd. Robust stability analysis is first investigated in this section, the outcome of which results in a feasible region λd ∈ [λ*d ,∞]such that the closed-loop system is stable. This provides the initial value of λd* for our automatic MPC tuning procedure. 3.1. General Structure. To aid the analysis, we have rearranged the closed-loop system in Figure 2. In this block diagram, the input uref and the disturbance d have been removed, as they do not affect the stability of the overall system; moreover, the nominal model has been combined with Fd for brevity, which is denoted as Fd = [F1d,F2d], where F1d and F2d are the transfer functions from u to yref and y to yref, respectively. In Figure 2, the real process is represented by P, which is different from the nominal model and can be written in the multiplicative uncertainty form P = G0(1 + W Δ)

M = −(1 + FMPC F̲ d1 + FMPC F̲ d2G0)−1FMPC F̲ d2G0W

(19)

where FMPC is the transfer function from Yref to u(k) and can be derived from the unconstrained MPC control law Δu(k) = f MPC(Uref, Yref). 3.2. Robust Stability Condition. For the M−Δ system, the robust stability has been well investigated in literature.27 This system is robustly stable if and only if det(1 − M(jω)Δ) ≠ 0, ∀ ω , ∀ |Δ| ≤ 1

(20)

For SISO systems, the above condition 20 is equivalent to |M(jω)| < 1, ∀ ω , ∀ |Δ| ≤ 1. ⇔ |W (jω)T (jω)| < 1, ∀ ω (21) −td0s

where T(s) = (1/(λdTp0s + 1))e is the sensitivity function of the closed-loop system. 3.3. Finding λ*d . Using the condition in eq 21, we may get a feasible region of λd such that robust stability is guaranteed. However, in this procedure, we should know how to choose the weighting function W so that the obtained region of λd is not too conservative. From the expression of P, we get Δ = W−1(P − G0)G−1 0 . As a part of the robust stability condition, we have |Δ| ≤ 1, which implies that

(16)

where Δ indicates the model uncertainty and W is a weighting function. Note that both P in Figure 2 and G in Figure 1 represent the unknown real process dynamics, but the difference is that G considers the effect of parametric uncertainty, whereas P represents the effect of unstructured uncertainty. By pulling out Δ, the closed-loop system can be represented in the general structure26 (see Figure 3a). The expression of N and the transfer function of the closed-loop system are as follows:

|W −1(P − G0)G0−1| ≤ 1

or

|W | ≥ |(P − G0)G0−1| (22)

⎡ N11 N12 ⎤ ⎥ N=⎢ ⎣ N21 N22 ⎦

(17)

F(N , Δ) = N22 + N21Δ(I − N11Δ)−1N12

(18)

This suggests us a way to construct the weighting function W using the maximum value of |(P − G0)G−1 0 |, which is known as the multiplicative error.28 In the parameter tuning framework introduced in Section 2, we assume the model uncertainty is reflected in the model E

DOI: 10.1021/ie5050583 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

operators to easily tell the consequences of choosing a certain combination of λ and λd. 4.1. Proposed Visualization Technique. Mathematically, the performance visualization problem is the combination of two sequences of optimization problems: Problem 1 For all t = 1, 2, ..., N and the process parameters Tp, td, and g constrained to their uncertainty regions, solve the following problems respectively. Here N is chosen to be a large positive integer to characterize the asymptotic behavior of the closed-loop system.

parameters. With this consideration, for a FOPDT model with the parametric uncertainty of the form g0 − Δg ≤ g ≤ g0 + Δg Tp0 − ΔTp ≤ Tp ≤ Tp0 + ΔTp td0 − Δtd ≤ td ≤ td0 + Δtd

the multiplicative error becomes (P − G0)G0−1 =

ge(td0− td)s(Tp0s + 1) − g0(Tps + 1) g0(Tps + 1)

(23)

max

y(t )

s.t.

td ∈ [ t ̲ d , td̅ ]; Tp ∈ [ T̲ p , Tp̅ ]; g ∈ [ g̲ , g ̅ ];

g , td , Tp

Then, the weighting function W can be constructed following a similar argument as that in Hu et al.28 Now, we determine the feasible region of λd; the approach is to find all λd’s such that |T(jω)| is less than |W(jω)| over all frequencies. From the expression of T(s), we have 1 |T (jω)| = (λdTp0ω)2 + 1 (24)

Xa(k + 1) = AaXa(k) + Ba Δu(k), y(k) = CaXa(k) + d(k), Δu(k) = fMPC (Uref , Yref ), k = 1, 2, ..., t ; Xa(0) = 0;

Thus, increasing λd will reduce |T(jω)|as well as the bandwidth of T(s); see Figure 4. This means that the system will become more

(25)

and min

y(t )

s.t.

td ∈ [ t ̲ d , td̅ ]; Tp ∈ [ T̲ p , Tp̅ ]; g ∈ [ g̲ , g ̅ ];

g , td , Tp

Xa(k + 1) = AaXa(k) + Ba Δu(k), y(k) = CaXa(k) + d(k), Δu(k) = fMPC (Uref , Yref ), k = 1, 2, ..., t ; Xa(0) = 0;

(26)

Note that model−plant mismatch is considered in the above optimization problems, as Δu(k) is calculated according to the nominal process parameters. The intuition of the two sequences of optimization problems is to characterize the upper and lower envelopes of the set of possible responses y(t) to Yref caused by the model mismatch. Several critical issues, however, exist in finding the exact solutions to the problem: (1) y(t) is not necessarily a convex function of the optimization parameters g, Tp, and td (in fact, g and Tp are explicitly expressed in Aa, whereas td affects Aa implicitly by controlling the size of the matrix); (2) the complexity of the dependence of y(t) on g, Tp, and td increases with the increase of t; (3) as performance visualization is only a small step in the overall iterative tuning procedure, the computation time allowed to solve Problem 1 is very limited. All these issues call for an efficient low-complexity heuristic solution to the performance visualization problem. To do this, we first rewrite problem into an equivalent form. Problem 2 For all t = 1, 2, ..., N and the process parameters Tp, td, and g constrained to their uncertainty regions, solve these respectively

Figure 4. Bode magnitude plot of T with respect to different λd.

robust as λd increases. As a result, to get the feasible region of λd, we only need to find the minimum λd (namely, λ*d) so that the robust stability condition is satisfied, which can be automatically calculated by using a bisection algorithm; for example, the feasible region of λd in Figure 4 is [2.5, ∞).

4. PERFORMANCE VISUALIZATION AND FAST IMPLEMENTATION As the performance indices considered in eq 15 are difficult to be analytically calculated, the structure of the tuning problem in eq 15 becomes ambiguous. To overcome this difficulty, we consider an alternative way of evaluating these indices, namely, the performance visualization technique to be introduced in this section. The objective is to graphically characterize the envelopes of the responses of a set of systems whose parameters lie within the user specified uncertainty intervals subject to a step set-point or load disturbance change, given the values of λ and λd. This not only helps to solve the tuning problem but also is of important help to the end users, as the visualization technique allows the

max

h(t , g , td , Tp , Uref , Yref , Xa(0), d)

s.t.

td ∈ [ t ̲ d , td̅ ]

g , td , Tp

Tp ∈ [ T̲ p , Tp̅ ] g ∈ [ g̲ , g ̅ ] F

(27) DOI: 10.1021/ie5050583 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

the Karush−Kuhn−Tucker (KKT) necessary conditions,29 the optimal solution to a constrained optimization problem can be obtained by enumerating all possible combinations of active constraints (namely, constraints in which the inequalities become equalities) and solving the resultant simplified optimization problems with equality constraints using Lagrange multiplier methods. On the basis of this observation, an efficient heuristic is to assume that the optimizer is achieved in the set of vertices of the parameter polytope, which constitutes a subset of the whole set of all possible active constraint combinations, and therefore, the problems can be solved by comparing the values of the objective function for only 23 vertices in the polytopic parameter space, which results in the following algorithm.

and min h(t , g , td , Tp , Uref , Yref , Xa(0), d)

g , td , Tp

s.t.

td ∈ [ t ̲ d , td̅ ] Tp ∈ [ T̲ p , Tp̅ ] g ∈ [ g̲ , g ̅ ]

(28)

Here, h(t, g, td, Tp,Uref,Yref,Xa(0),d) = y(t) is used to explicitly express the dependence of y(t) on td, Tp, Uref, Yref, Xa(0), and d, and it is obtained by composing y(t) with the state space equation and the MPC law fMPC(Uref,Yref). In Problem 2, we observe that both optimization problems are nonlinear optimization problems within polytopes. According to

Remark 1 Although this heuristic is not guaranteed to be optimal, it is intuitive in that extreme behavior of the step responses mostly happens at the extreme process parameters. Also, this method is applicable to the characterization of envelopes for other process variables, for example, control variables. In addition, as the envelopes are obtained from the corresponding responses of the extreme processes, the algorithm can be applied to the cases of both MPC with more complex process structures and MPC with constraints; the cost is that the computational burden will be heavily increased, as a real time MPC is needed to simulate the process responses. This further interprets why an unconstrained MPC is considered, for which the control gain can be analytically calculated off-line to speed up the algorithm. Note that another important reason to choose the unconstrained MPC specifically in this work is that the constraints on the inputs are relatively loose for most operating conditions of industrial paper machines and, thus, can be ignored so as to solely focus on the difficulties introduced by directly considering the user-specified parametric uncertainty and the performance specifications on overshoots and settling times. 4.2. Numerical Evaluation of Time-Domain Performance Indices. One of the major benefits of the above proposed performance visualization technique is that for arbitrary userspecified parameter uncertainty intervals, it allows direct evaluation of time-domain performance indices, the expressions of which are normally not possible to calculate analytically. This provides the feedback information for the overall iterative parameter tuning procedure. Assuming all responses in the set have the same final values [because the implemented process in eq 4 can track step setpoint changes, this assumption is automatically satisfied], we introduce the definitions of worst-case overshoot and settling time for a set of step responses, which are generalized from the definitions of overshoot and settling time for a single response. Definition 1 (Worst-case overshoot) The worst-case overshoot ζw of a set of step responses with the same final value is the maximum value in all responses minus the final value divided by the final value.

Definition 2 (Worst-case settling time) The worst-case settling time ts of a set of step responses with the same final value is the time required for all the responses to reach and stay within a range of prespecified percentage [the 5% criterion is assumed throughout this work] of the final value. Based on these definitions and the proposed visualization method, the values of ζw and ts can be calculated numerically. The worst-case overshoot ζw can be computed by finding the maximum peak in y(t) and max{y ̅ (t )} − y(∞) ζw =

t

y(∞)

(29)

while the numerical evaluation of the worst-case settling time can be implemented by reversing the time order of sequences y and y, and identifying the time tag of the first element that escapes the ±5% interval. 4.3. Numerical Verification of the Visualization Procedure. In this subsection, we illustrate the developed visualization technique through extensive simulations. The proposed procedure is applied to typical FOPDT processes in process control, and the results are shown in Figure 5. Three types of processes are considered, with the first one having balanced time constant and delay (see Figure 5a−c), the second one being delay dominant (see Figure 5d−f), and the third one being time-constant dominant (see Figure 5g−i). For each type of process, the values of (λ, λd) are set to (1.5, 1.5), (2, 1), and (1, 2), respectively. To demonstrate the efficiency of the visualization procedure, the responses of 100 randomly generated systems satisfying the parameter uncertainty intervals are also plotted in each case. Due to the space limitation, the parameters of each case and the measured overshoots and settling times are listed in the corresponding subfigures of Figure 5. From these figures, it is observed that the resultant envelopes can always efficiently characterize the set of possible responses generated by the processes with parameters lying within the specified intervals, even when the systems are unstable (see Figure 5e); the numerical evaluation of settling times and overshoots is shown to be accurate. Also, notice G

DOI: 10.1021/ie5050583 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 5. Numerical verification of the proposed performance visualization technique.

5. ITERATIVE PARAMETER TUNING In this section, an iterative tuning procedure is proposed for λ and λd, based on the stabilizing region of λd determined in Section 3.3 and the proposed visualization technique in

that the computation complexity of the procedure is satisfactory, as a single run of the procedure for a FOPDT process takes only 0.17 s on a laptop with Intel core-i5 and 6 GB memory. H

DOI: 10.1021/ie5050583 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 6. Numerical verification of the properties of ts(λ, λd) and ζw(λ, λd). To improve visibility, large values of overshoot and settling time are truncated without affecting the verification results, and instead of the value of overshoot, the logarithm of overshoot is plotted to verify the monotonicity property of ζw(λ, λd).

In addition, given the stability of the system (i.e., λd > λ*d), the worst-case overshoot ζw(λ, λd) will empirically decrease toward 0 as λ → ∞, following from the above monotonic property of ζw(λ, λd) with respect to λ. This implies that as long as ζ*w > 0, there always exists a pair (λ, λd) such that the constraint in the tuning eq 15 is satisfied, and thus the feasibility of the tuning problem can be ensured in general. Also, note that in the above feasibility analysis, the effect of parametric uncertainty level is reflected in the value of λ*d; however, this does not affect the feasibility of the tuning problem, as there is no restriction on the feasible region of λ and λd. 5.2. Line-Search Based Iterative Tuning Algorithm. Based on the empirical properties of ts(λ, λd) and ζw(λ, λd), an automatic tuning algorithm is proposed in this subsection. For industrial MPC tuning, the requirements on the efficiency and the reliability of the algorithms are emphasized, because (1) a satisfactory solution obtained within a small amount of time results in better user experience, compared with that of an optimal solution obtained for certain models (with unavoidable modeling errors) at the cost of more computation time; (2) the explicit expressions of ts(λ, λd) and ζw(λ, λd) are unknown and only numerical values are available using the visualization technique, which limits the allowed type of algorithms that can be considered here, for example, Newton and quasi-Newton algorithms are prohibitive, because the numerical evaluation of the first order derivatives could lead to unexpected errors and the failure of the overall tuning procedure. Based on these considerations, we choose line-search methods29 in finding a satisfactory combination of λ and λd. For the proposed algorithm (see Algorithm 2), the nominal process parameters Tp0, td0, and g0 are assumed to be known, which can be readily obtained by industrial process identification software, for example, the Honeywell Profit Design Studio; the uncertainty on the process parameters are assumed to be known as well, which are specified by the end users. The algorithm is performed iteratively to find a pair of λ and λd; each iteration is composed of three stages, which are explained in detail in the following. The algorithm stops either when a stationary point is arrived (i.e., the same pair of λ and λd is obtained in two consecutive iterations) or when the computation time runs out (which is counted as the number of iterations allowed), and it collects the pair of feasible (λ, λd) with the smallest settling time during the iterative procedure.

Section 4.1. For notational simplicity, we explicitly use ts(λ, λd) and ζw(λ, λd) to represent the relationship of worstcase settling time and worst-case overshoot with λ and λd, respectively. 5.1. Empirical Unimodality/Monotonicity of ts(λ, λd) and ζw(λ, λd). In order to develop algorithms to find a pair (λ, λd) that solves the problem in 15, geometric properties of ts(λ, λd) and ζw(λ, λd) have to be explored. As the analytical expressions of ts(λ, λd) and ζw(λ, λd) in general do not exist, the remedy taken here is to first provide a qualitative analysis on the properties of the functions and then to numerically verify these properties. According to Section 3, the stability of the system is determined by λd and the robustness of the system increases with the increase of λd. Therefore, it is intuitive from an engineering perspective that ts(λ, λd) is a unimodal function of λd for a fixed value of λ: A very small value of λd causes a large settling time, due to the relative aggressive and oscillatory response; a large value λd causes a large settling time as well, due to the over sluggish responses. Similarly, when λd is fixed, ts(λ, λd) is also a unimodal function of λ. On the other hand, the Fr filter controlled by λ controls only the speed of the response, and thus, a larger value of λ leads to a smaller value of overshoot. In this way, ζw(λ, λd) can be empirically treated as a monotonically decreasing function of λ. To verify this analysis, we numerically evaluate the relationship of ts(λ, λd) and ζw(λ, λd) with λ and λd for different FOPDT processes. The envelopes in the visualization procedure are generated within Tp ∈ [(1 − r)Tp0, (1 + r)Tp0], td ∈ [(1 − r)td0, (1 + r)td0], g ∈ [(1 − r)g0, (1 + r)g0], r being the uncertainty level. Again, three types of typical processes are considered, with the first one having balanced time constant and delay, the second one being delay dominant, and the third one being time-constant dominant. For each type of processes, the levels of parameter uncertainty are set to r = 10%, r = 30%, and r = 60%, respectively, to consider different uncertainty level combinations. The plots of the relationship of ts(λ, λd) and ζw(λ, λd) with respect to λ and λd for Tp0 = 90 s, td0 = 90 s, g0 = 0.129, and r = 30% are shown in Figure 6, where the empirical properties of ζw(λ, λd) and ts(λ, λd) can be observed. As the plots of ζw(λ, λd) and ts(λ, λd) for other parameter settings are observed to have similar shapes and properties, they are safely omitted here considering the space limitation. I

DOI: 10.1021/ie5050583 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

overshoot with respect to λ, this can be achieved by utilizing a bisection algorithm, which leads to the codes in lines 22−30. In the third stage (lines 31−42), the settling time is further optimized with respect to λ within the region λ ∈ [λ,∞), λ being calculated in the second stage and λd being calculated in the first stage. Notice that due to the monotonicity property of ζw(λ, λd) with respect to λ, the overshoot constraint in eq 15 is satisfied for all λ ≥ λ. Within this region, ts(λ, λd) is either a unimodal function or a monotonic function of λ, which allows the golden section search algorithms to find the λ that achieves the smallest settling time. If the stopping conditions are not satisfied, the algorithm proceeds to stage 1, with the value of λ setting to that obtained at this stage (i.e., stage 3).

In the first stage (lines 10−21), λd is first tuned to get the optimal settling time for a fixed λ. At this stage, λ is set to a very small value, which would lead to a relatively aggressive response for a fixed λd. By the empirical unimodality property, the optimization of the settling time with respect to λd is achieved by golden section search with almost linear convergence rate [guaranteed linear convergence rate can be achieved as well by using Fibonacci sequences, as well29] without requiring the numerical calculation of the derivatives. In the second stage (lines 22−30), λ is further tuned to find λ that activates the constraint on overshoot (namely, the constraint becomes an equality constraint), based on the tuned value of λd in stage 1. From the empirical monotonicity property of the J

DOI: 10.1021/ie5050583 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 7. Performance trade-off for uncertainty level [−50%, 100%].

Figure 8. Performance trade-off for uncertainty level [−50%, 50%].

⎛⎡ ⎤ ⎡ ⎤ ϵ ϵ ⎜⎜⎢log ⎥ + ⎢log 0.5 ⎥⎥ 0.618 ⎢ * ̲ 100 0.02 − 100 − λ ⎥ ⎢ ⎝ d

Remark 2 Since line search methods are used in the three stages of the algorithm, the computational complexity can be approximately evaluated. Following from the principle of line search, the number of inner iterations in stage 1 is bounded by ⌈ log0. 618(ϵ/(100−λ*d))⌉, where logab denotes the logarithm of b with base a, and ⌈ x ⌉ denotes the smallest integer that is greater than or equal to x. Similarly, the numbers of inner iterations in stage 2 and stage 3 are bounded by⌈ log0. 5(ϵ/ (100−0.02))⌉ and ⌈ log0. 618(ϵ/100)⌉, respectively. Let τv denote the time required to calculate the worst-case overshoot and settling time for a pair (λ, λd). Noticing that the same operation, namely, calculating the worst-case overshoot and settling time for a given pair (λ, λd), is performed in each inner iteration of the three stages and that the number of outer iterations is bounded by N0 (see line 7 of Algorithm 2), the total computation time of the algorithm is upper bounded by

⎡ ϵ ⎤⎞ + ⎢log 0.618 ⎥⎟⎟N0τv ⎢ 100 ⎥⎠

Remark 3 If the constraints in the MPC formulation are not ignored for parameter tuning, the empirical monotonicity and unimodality properties still hold in general, and thus Algorithms 1 and 2 and the discussions in Remark 2 remain valid. The major difference, however, is that a real-time MPC solver has to be used in Algorithm 1 to generate the responses yi(t) for t = 1, 2, ..., N (typical values of N lie between 100 and 400) and i = 1, 2, ..., 8, which significantly increases τv, namely, the time needed to calculate the worst-case overshoot and settling time. To see this, note that for the unconstrained case, τv is approximately 0.2 seconds on a laptop with Intel core-i5 and 6 GB memory, while K

DOI: 10.1021/ie5050583 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 9. Performance trade-off for uncertainty level [−25%, 30%].

Table 2. Tuning results for uncertainty level [−50%,100%] overshoot specification ζ*w proposed algorithm

brutal search

ζw ts λ λd ζw ts λ λd

10%

20%

30%

40%

50%

10% 2505 s 8.8279 3.61 9.1214% 2535 s 9.0214 3.62

16.98% 2415 s 8.0097 3.7413 16.98% 2415 s 8.0097 3.7413

29.90% 2340 s 6.6243 3.9537 29.90% 2340 s 6.6243 3.9537

33.42% 2340 s 6.2889 4.0039 30.3078% 2340 s 6.5981 4.095

44.6159% 2355 s 5.3149 4.1662 44.6159% 2355 s 5.3149 4.1662

predicting the process response yi(t) for a single step t using an industrial MPC solver on the same computer takes up 0.03 seconds.

specifications on worst-case overshoots. To test the performance of the algorithm for optimality in worst-case settling times, brutal search is performed for each level of uncertainty and each overshoot specification by enumerating all points contained in a prespecified gridded region of the λ parameters, which finds the almost-optimal solution to the optimization problem in eq 15 with closeness to optimality controlled by the resolution of the grid. The obtained performance comparison is presented in Figure 7−9, where the computational time of the proposed algorithm is also indicated for each point. It is observed that although the outcome of the proposed algorithm does not have theoretically guaranteed optimality, it is always close to the result of the brutal search, which normally takes 10 min to calculate the best pair of (λ, λd) for a single specification of worst-case overshoot. On the other hand, the computation time of the tuning algorithm is also satisfactory, which normally takes 3−5 s, depending on the size of uncertainty. To take a closer look on the tuning parameters and performance indices, the tuning results for uncertainty level [−50%, 100%]and overshoot specifications 10%, 20%, ..., 50% are further presented in Table 2. In addition, the obtained closed-loop step response envelope for uncertainty level [−25%, 30%]and ζ*w = 20% is shown in Figure 10, and the corresponding envelope for control signals is shown in Figure 11 to further illustrate the tuning results. Now we test the tuning design results in the Honeywell real time MPC + Simulator environment. To account for model mismatch, the real process is taken as

6. INDUSTRIAL EXAMPLE In this section, we illustrate the effectiveness of the proposed tuning algorithm with the results obtained for a typical MD process in a paper machine. Consider the following nominal model that is used in machine directional paper machine control: G0(s) =

0.0135 −90s e 60s + 1

(30)

which is identified by industrial process identification software based on input−output data from an industrial site and models the dynamics from stock to conditioned weight. The system is discretized at sampling time Ts = 15 s. The MPC weighting matrices are set to Q1 = I, Q2 = 0, and Q3 = 0, the prediction horizon Hp is 42 and the control horizon Hu is 20. For notational brevity in this section, we represent the uncertainty level in a compact form [−r %, r ̅ %], if the parametric uncertainties on the nominal process parameters are specified as Tp ∈ [(1 − r̲ %)Tp0 , (1 + r ̅ %)Tp0] td ∈ [(1 − r̲ %)td0 , (1 + r ̅ %)td0] g ∈ [(1 − r̲ %)g0 , (1 + r ̅ %)g0]

For the nominal model in eq 30, the proposed tuning algorithm is applied for different levels of parameter uncertainty and different L

DOI: 10.1021/ie5050583 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 10. Step response envelope for uncertainty level [−25%, 30%]and ζ*w = 20%.

Figure 11. Control input envelope for uncertainty level [−25%, 30%]and ζ*w = 20%.

G (s ) =

0.0174 −116.1s e 77.4s + 1

overshoot is taken as 20%. Three sets of λ parameters obtained by the proposed tuning algorithm for different uncertainty levels, namely, [−25%, 30%], [−50%, 50%] and [−50%, 100%], are applied to the process simulator and Honeywell real-time MPC for performance comparison. The obtained responses are shown in Figure 12. It is observed that despite the model− plant mismatch and measurement noise, the system outputs for all the three sets of λ parameters always robustly track the setpoint for all changes of the operating conditions, which further validates the effectiveness of the MPC tuning parameters obtained from the proposed algorithm. Also, it is observed that the λ parameters obtained for higher uncertainty levels (e.g., those for uncertainty level [−50%, 50%]) lead to relatively smoother responses, compared with that for uncertainty level [−25%, 30%]. The interpretation is that if the parameters are designed to tolerate more uncertainties, the corresponding control actions have to be less aggressive to maintain performance (in terms of overshoots) for all possible processes lying within the uncertainty set around the nominal process, and thus, the responses are normally smoother compared with

(31)

which lies within the uncertainty level [−25%, 30%] of the nominal process in eq 30. The initial operating point is y(0) = 432, u(0) = 3790, according to the actual operating conditions. The MPC optimization parameters Hp, Hu, Q1, Q2, and Q3 are chosen to be the same as those used for the above tuning results, and the following constraints on the control signals are incorporated in MPC implementation 0.9u(0) ≤ U ≤ 1.1u(0) −0.1u(0) ≤ ΔU ≤ 0.1u(0)

To consider possible changes of operating conditions, a set-point change of 2 lbs/1000 ft2 is made at t = 300 s, an output disturbance of 2 lbs/1000 ft2 is in effect at t = 2250 s, an input disturbance of 80 gpm is introduced at t = 4500 s, and the measurement noise is taken to be zero-mean Gaussian with standard deviation 0.1 lbs/1000 ft2. The specification on worst-case M

DOI: 10.1021/ie5050583 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 12. Real time MPC + Simulator results.



those obtained for optimal tracking, which are less robust to uncertainties.

(1) Shi, D., Wang, J., Forbes, M., Backström, J., Chen, T. Method and apparatus for specifying and visualizing robust tuning of model-based controllers. U.S. Patent Application No. 61/954.912, 2014. (2) Backström, J., Baker, P. A Benefit Analysis of Model Predictive Machine Directional Control of Paper Machines Pulp Pap. Can. 2008. (3) Chu, D., Forbes, M., Backström, J., Gheorghe, C., Chu, S. Model Predictive Control and Optimization for Papermaking Processes. In Advanced Model Predictive Control; Zheng, T. Ed.; InTech: Rijeka, Croatia, 2011. (4) Åström, K.; Hägglund, T. The future of PID control. Control Eng. Pract. 2001, 9, 1163−1175. (5) Mayne, D.; Rawlings, J.; Rao, C.; Scokaert, P. Constrained model predictive control: Stability and optimality. Automatica 2000, 36, 789− 814. (6) Rawlings, J., Mayne, D. Model Predictive Control: Theory and Design, 1st ed.; Nob Hill Publishing: Madison, WI, 2009. (7) Korda, M.; Gondhalekar, R.; Oldewurtel, F.; Jones, C. Stochastic MPC Framework for Controlling the Average Constraint Violation. IEEE Trans. Autom. Control 2014, 59, 1706−1721. (8) Narang, A.; Shah, S.; Chen, T.; Shukeir, E.; Kadali, R. Design of a model predictive controller for interface level regulation in oil sands separation cells. Am. Control Conf. 2012, 2812−2817 2012. (9) Jiang, H.; Shah, S.; Huang, B.; Wilson, B.; Patwardhan, R.; Szeto, F. Model analysis and performance analysis of of industrial MPCs. Control Eng. Pract. 2012, 20, 219−235. (10) Zhou, K., Doyle, J., Glover, K. Robust and Optimal Control. Prentice Hall: Upper Saddle River, NJ1996. (11) Vlassis, N.; Jungers, R. Polytopic uncertainty for linear systems: New and old complexity results. Syst. Control Lett. 2014, 67, 9−13. (12) Al-Ghazzawi, A.; Ali, E.; Nouh, A.; Zafiriou, E. On-line tuning strategy for model predictive controllers. J. Process Control 2001, 11, 265−284. (13) Garriga, J.; Soroush, M. Model predictive controller tuning via eigenvalue placement. Am. Control Conf. 2008, 429−434. (14) Di Cairano, S.; Bemporad, A. Model Predictive Control Tuning by Controller Matching. IEEE Trans. Autom. Control 2010, 55, 185− 190. (15) Shah, G.; Engell, S.; Tuning MPC for desired closed-loop performance for SISO systems. 2010 18th Mediterranean Conference on Control & Automation (MED); IEEE: Piscataway, NJ, 2010; pp 628− 633. (16) Shah, G.; Engell, S. Tuning MPC for desired closed-loop performance for MIMO systems. Am. Control Conf. 2011, 4404−4409.

7. CONCLUSION In this work, a systematic procedure is introduced for automatic parameter tuning of 2-DOF MD-MPC for time-domain robust performance under parametric uncertainties. An efficient technique for time-domain closed-loop performance visualization in the presence of parametric model−plant mismatch specifications has been proposed. Under these model−plant mismatch specifications, an algorithm to automatically determine the controller tuning parameters λ and λd that exactly satisfy the specification on worst-case overshoot while guaranteeing a satisfactory worstcase settling time is proposed. The method has been designed for quality control of paper machines but is equally applicable to general industrial processes. Future research topics include the consideration of variations in process outputs as an additional performance index, as well as techniques to further increase the speed of the tuning algorithm.



REFERENCES

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Present Address

§ State Key Laboratory of Intelligent Control and Decision of Complex Systems, School of Automation, Beijing Institute of Technology, Beijing 100081, China.

Notes

The authors declare no competing financial interest. A portion of the results presented in this work was submitted as a provisional U.S. patent.1



ACKNOWLEDGMENTS This work was supported by NSERC and Honeywell. The authors thank Ning He from the University of Alberta for helping implement part of the codes of this work. The suggestions of Dr. Greg Stewart from Honeywell Automotive Software are also gratefully appreciated. Finally, the authors would like to thank the Associate Editor and the anonymous reviewers for their suggestions, which have improved the quality of the work. N

DOI: 10.1021/ie5050583 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research (17) Toro, R.; Ocampo-Martínez, C.; Logist, F.; Impe, J. V.; Puig, V. Tuning of Predictive Controllers for Drinking Water Networked Systems. 18th IFAC World Congr. 2011, 14507−14512. (18) Vallerio, M.; Impe, J. V.; Logist, F. Tuning of NMPC controllers via multi-objective optimization. Comput. Chem. Eng. 2014, 61, 38−50. (19) Garriga, J. L.; Soroush, M. Model Predictive Control Tuning Methods: A Review. Ind. Eng. Chem. Res. 2010, 49, 3505−3515. (20) Bordons, C.; Camacho, E. A generalized predictive controller for a wide class of industrial processes. IEEE Trans. Control Syst. Technol. 1998, 6, 372−387. (21) Kong, H.; Goodwin, G.; Seron, M. Predictive metamorphic control. Automatica 2013, 49, 3670−3676. (22) Chu, D., Forbes, M., Backström, J. Technique for converting a model predictive control (MPC) system into an explicit two-degrees of freedom (2DOF) control system. U.S. Patent Application No. 13/ 907495, 2013. (23) Fan, J. Model Predictive Control for Multiple Cross-Directional Processes: Analysis, Tuning and Implementation. Ph.D. Thesis, The University of British Columbia, 2003. (24) Fan, J.; Stewart, G.; Dumont, G.; Backström, J.; He, P. Approximate steady-state performance prediction of large-scale constrained model predictive control systems. IEEE Trans. Control Syst. Technol. 2005, 13, 884−895. (25) Hang, C.; Chin, D. Reduced order process modelling in selftuning control. Automatica 1991, 27, 529−534. (26) Zhou, K., Doyle, J. Essentials of robust control; Prentice-Hall: Englewood Cliffs, NJ, 1998. (27) Skogestad, S., Postlethwaite, I. Multivariable Feedback Control: Analysis and Design; John Wiley & Sons: New York, 1996. (28) Hu, J.; Bohn, C.; Wu, H. Systematic H∞ weighting function selection and its application to the real-time control of a vertical take-off aircraft. Control Eng. Pract. 2000, 8, 241−252. (29) Bazaraa, M., Sherali, H., Shetty, C. Nonlinear Programming: Theory and Algorithms, 3rd ed.; John Wiley & Sons, Inc.: Hoboken, NJ, 2006.

O

DOI: 10.1021/ie5050583 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX