An Algorithm for Tuning NMPC Controllers with Application to

Subscriber access provided by United Arab Emirates University | Libraries Deanship

Article

An Algorithm for Tuning NMPC Controllers with Application to Chemical Processes. Federico Lozano Santamaría, and Jorge M. Gómez Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.6b01121 • Publication Date (Web): 05 Aug 2016 Downloaded from http://pubs.acs.org on August 10, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

An Algorithm for Tuning NMPC Controllers with Application to Chemical Processes. Federico Lozano Santamaría and Jorge M. Gómez*. Grupo de Diseño de Productos y Procesos, Departamento de Ingeniería Química, Universidad de los Andes. Carrera 1 No. 18A-10, Bogotá, Colombia. Key words NPMC tuning, economic oriented NMPC, NLP sensitivity, utopia tracking. Abstract The advantages of linear/nonlinear model predictive control (N)MPC for dealing with the multiple input multiple output problem, for performing optimization and for handling constraints are well known and because of that it has been applied widely in the chemical industry. However, there is a recurrent problem for this kind of controllers and it is how to define the best tuning parameters to achieve a good closed-loop response. This is an open question for research even when the common practice is to define the (N)MPC tuning parameters by trial and error or by the expertise of the control engineer. There are some works that propose a systematic

*

To whom correspondence should be addressed. E-mail: [email protected]

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 45

definition of the (N)MPC tuning parameters, but most of them are restricted to certain prediction models or do not define all the tuning parameters. This papers presents a general tuning algorithm for (N)MPC controllers that allows the systematic definition of all the tuning parameter and it is not restricted to a specific prediction model or to a problem formulation. It is based on multiobjective optimization for the definition of the objective function weights and on the optimization of a closed loop performance index for the definition of the horizon lengths. Three cases of chemical processes are considered in detail to illustrate the application and advantages of the proposed (N)MPC tuning algorithm.

Nomenclature Subscript 0 : initial condition. Inlet feed condition for case 1 and case 2. 𝐷 : distillate 𝑓 : variable evaluated at the final time of the prediction horizon. 𝑘 : index used to define a discrete point. 𝑝 : variable associated with an artificial constraint for the parameter 𝑝 𝑤 : coolant fluid Superscript 𝐸𝑡𝑂𝐻 : ethanol 𝐿 : lower bound 𝑠𝑝 : set-point 𝑈 : upper bound. Latin symbols 𝑎 : parameters of the response surface that describe a KPI. 𝐴 : transpose Jacobian of 𝑐 𝑐: vector of equality constraints in an optimization problem


2

Page 3 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


𝐶 : a constant. Molar concentrations in case 1 and case 2. ̅̅̅ 𝐶𝐶 : average molar concentration of propylene glycol (C) in case 1 and case 2. 𝐶𝐷 : distillate profit in case 3. 𝐶𝑝 : specific heat capacity. 𝐶𝑄 : energy cost in case 3. 𝐸 : activation energy of a reaction rate. 𝐹 : molar flow rate. 𝐻𝐿 : Hessian matrix of the Lagrange function 𝐼: identity matrix. 𝑘0 : frequency factor of a reaction rate. 𝐾 : key performance index function ̂𝑖 : key performance index value 𝐾 𝐿 : Lagrange function. 𝑚̇ : mass flow rate. 𝑁𝑝 : prediction horizon. 𝑁𝑢 : control horizon. 𝑝 : parameter of an optimization problem. 𝑄 : heat transferred. 𝑄𝑅 : reboiler duty. 𝑟 : reaction rate. 𝑅 : gas constant in case 1 and case 2. Molar reflux ratio in case 3. 𝑡 : time 𝑡𝑠𝑜𝑙 : maximum time for solving the OCP problem 𝑇 : temperature 𝑢 : manipulated variable. 𝑈𝐴 : global heat transfer coefficient times the surface area. 𝑣: vector of dual variables in a optimization problem.


3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 45

𝑉 : diagonal matrix with the values of 𝑣. Reactor volume for case 1 and case 2. 𝑤𝐼 : disturbance in the input states 𝑤𝑂 : disturbance in the output states 𝑊 : Hessian matrix of the Lagrange function. 𝑊𝑖 : weight for a specific objective in the OCP objective function. 𝑥 : states. Molar fraction in case 3. 𝑥𝑎 : algebraic states 𝑥𝑑 : differential states ̅̅̅̅̅̅̅̅ 𝑥𝐷𝐸𝑡𝑂𝐻 : average ethanol molar percentage in distillate for case 3. 𝑦 : control outputs. 𝑧 : vector of all the variables (𝑥, 𝑦) in a nonlinear optimization problem. 𝑍 : diagonal matrix with the values of 𝑧. Greek symbols Δ𝐻𝑟𝑥𝑛 : reaction enthalpy. 𝝓 : vector of objective functions. 𝜙𝑖 : a particular objective 𝑖 / without the index is an objective function 𝜆 : Lagrange multiplier for the equality constraints. 𝜆̅ : Lagrange multiplier for the artificial constraints. 𝜏 : sampling time. 𝜗 : estequiometric coefficient.


4

Page 5 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction Nonlinear model predictive control (NMPC) is an advanced control strategy widely used in different industries because of its capability to handle constraints, to perform optimization and to handle the multivariable problem directly. The essence of NMPC is to solve an optimal control problem (OCP) in a receding horizon scheme to construct the control law, in other words, at every time period it minimizes certain performance index subject to a (non)linear model (equality constraints) and operational limits (inequality constraints). Even though (N)MPC controllers have been used in many applications it is clear that their good performance depends on two important factors: the use of an accurate enough model to represent the plant dynamic and the selection of tuning parameters to obtain a good closed loop response 1. If the controller exhibits poor performance, even when its tuning parameters have been improved, it is possible that the model used for prediction is not appropriate and it should be refined 1. This paper is not focused on the selection of prediction models for (N)MPC, rather it assumes that a nonlinear model can describe the process dynamic precisely. On the other hand, if the model is adequate, the closed loop performance of the system can be improved selecting appropriate tuning parameters for the controller. Usually, the (N)MPC tuning has been done via trial and error or based on the designer expertise 2, but it is possible to tune the controller systematically and there are many studies that try to solve the tuning problem from different perspectives 1. In addition, the increasing number of (N)MPC applications stipulates a systematic tuning procedure. The traditional trial and error methodology for tuning (N)MPC controllers is an inefficient strategy because it requires the evaluation of a large number of combinations. Also, many of the parameters have coupled or overlapping effects on the performance of the controller; thus, it is clear the need for a more precise and consistent method for tuning these controllers. The main


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 45

idea of a systematic tuning procedure for (N)MPC controllers is to define the controller parameters (e.g. prediction horizon, control horizon, weights of the different objectives) based on performance indexes of the closed loop behavior of the system. It might seem easy to define a performance index for the closed loop response as a similar index is defined for the OCP that is solved periodically, but there could be more aspects that are important for the control engineer that are not represented in the objective function of the OCP. Moreover, the OCP is usually nonlinear and/or is subject to inequality constraints which make it impossible to find an analytical expression of the control law 3. In consequence, the performance index defined for the closed loop behavior of the system cannot be related explicitly with the tuning parameters of the controller. In the literature there are many works dealing with the problem of tuning (N)MPC controllers from different perspectives, but none of them, as far as we know, is capable to define all the tuning parameters for any kind of (N)MPC controller and the development of a more general methodology is still a challenge and an open area for research 4. For example, Shridhar and Cooper 5 present a tuning strategy based on a compilation of heuristics and on the definition of a well-conditioned optimization problem. This strategy begins with the approximation of the prediction model to a first order plus dead time (FOPDT) model and then it uses some simple equations to calculate all the tuning parameters of the controller, including the length of the horizons. Nevertheless, the approximation to a FOPDT model and the use of heuristics might not be accurate enough to guarantee a good closed loop performance or a robust controller. It is important to note that this approach can also be applied to nonlinear MPC, but first the model must be approximate to a FOPDT model in order to identify some characteristics of the process such as time constants and dead times which are required in the tuning expressions. A more


6

Page 7 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


rigorous approach proposed by Di Carino and Bemporad 6 is a tuning strategy based on matching the MPC closed loop performance to the closed loop performance of a “favorite” controller. In order to match the two controllers they propose a minimization of the difference of the control law predicted by both controllers. However, it only defines the objective function weights and it is limited to linear models without inequality constraints in which it is possible to find an explicit solution of the feedback control sequence. In another work by Francisco and Vega

7

a

multiobjective optimization strategy, specifically goal programming, is used to solve the conflicting objectives of the controller and to define the penalization weights of the rate of change of the manipulated variables. Then, a random search is used to define the length of the horizons based on the same performance index of the multiobjective optimization problem. This strategy is limited to linear models without constraints and for nonlinear models it is necessary to linearize them around the operation point which reduces the efficacy of the final tuning parameters on the closed-loop performance. One final example is the tuning procedure presented by Vallerio et. al.

8

in which a multiobjective optimization procedure is used to define the

weights of the conflicting objectives assuming a fix control and prediction horizon. The multiobjective OCP is solved using the enhanced normalized normal constraint or normal boundary intersection method and the Pareto front is defined. Once the analyst has chosen a compromised solution on the Pareto front, the result of the multiobjective optimization can be translated into a weighted sum which is a simpler and more common approach for NMPC controllers. In this paper we present a tuning algorithm for NMPC controllers that allows the definition of the tuning parameters (prediction horizon, control horizon and the objective function weights) and that is not restricted to a specific type of prediction model or problem formulation. The


7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 45

proposed tuning algorithm is based on the constraint optimization of a performance index of the closed loop behavior of the controller. It uses a multiobjective optimization technique, specifically a utopia point tracking approach, to define the weights of the different objectives and a nonlinear fit to approximate the behavior of the closed-loop performance indexes with respect to the length of the horizons. This final step allows the formulation and solution of a new optimization problem in which the tuning parameters of the controller are the optimization variables. The rest of this paper is structured as follows: Section 2 presents a brief description of NMPC and dynamic optimization. The tuning strategy, which is the main contribution of this work, is described in detail in Section 3. It presents each step of the algorithm and the different optimization techniques used within it such as multiobjective optimization using the utopia tracking approach and the approximation of the performance indexes to quadratic functions of the prediction and control horizons. Section 4, describes the advanced step NMPC (asNMPC) controller. Section 5, illustrates the application of the NMPC tuning algorithm with two case studies of chemical processes, a CSTR reactor and an extractive distillation column. For the first case, two possible scenarios are considered which result in a final repertory of three examples that apply the proposed tuning algorithm. The final conclusions and considerations of the paper are presented in Section 6.

2. Background on NMPC and dynamic optimization Nonlinear model predictive control (NMPC) is an advanced control strategy that solves an OCP at every time period. The solution to this problem is a control sequence over a future horizon, but from this sequence only the first element is implemented in the plant to construct the control law.


8

Page 9 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1 presents a simple example of how model predictive control works. At the current time (𝑡𝑘 ) an OCP subject to a model of the process is defined and solved online in a future prediction horizon (𝑁𝑝 ) where only a number of control actions are allowed (𝑁𝑢 ). The solution of the optimization problem is the prediction of the future behavior of the process and the necessary control actions to achieve the desired operation. However, only the first value of the manipulated variables is feedback to the process and the rest future actions predicted are discarded. The definition of the OCP is updated based on new measurements of the state variables and it is solved again. Past moves and trajectory

Current time

Past output measurements Predicted optimal trajectory Set-point Current state Optimal control sequence Feedback value

Control horizon [k+Nu]

Prediction horizon [k+Np]

Past control actions (Manipulated variable history) t(k-1)

t(k)

t(k+1)

t(k+2)

Time

Figure 1. Representation of the receding horizon strategy in the (N)MPC formulation The critical step in the implementation of NMPC controllers is the solution of the OCP because it has to be solved periodically and in less time than a sampling time interval. A general formulation of the OCP is given by Eq. (1), where 𝑥𝑑 are the differential states, 𝑥𝑎 the algebraic states, 𝑦 the measured outputs, 𝑤𝐼 the disturbances in the states, 𝑤𝑂 the disturbances in the outputs and 𝑢 the manipulated variables. The objective function Eq. (1a) is composed by two terms: a running cost (𝜙) which is a weighted sum of different objectives like a quadratic


9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 45

tracking function, a quadratic penalization of the rate of change of the manipulated variables and/or an energy or economic function; and a terminal cost (𝜙𝑓 ) that is only a function of the states at the final time of the prediction horizon and is normally used to approximate the finite horizon problem to the infinite horizon problem 9. The optimization problem is subject to a set of differential and algebraic equations Eqs. (1b) – (1d) that can be linear defining a quadratic programming problem (QP) or nonlinear defining a nonlinear programming problem (NLP). The Eq. (1f) is also an equality constraint, but it is presented separately from the others because it represents the measured outputs for the control problem. The inequalities of the problem, Eq. (1e), are normally associated with variables bounds or bounds on their rate of change, but they can also be nonlinear expressions. 𝑡𝑘 +𝜏𝑁𝑝

min 𝑂𝐵𝐽 = ∫

𝑥𝑑 ,𝑥𝑎 ,𝑢

(1a)

𝜙(𝑡, 𝑦, 𝑢)𝑑𝑡 + 𝜙𝑓 (𝑦𝑓 , 𝑢𝑓 )

𝑡𝑘

𝑑𝑥𝑑 𝑠. 𝑡. = 𝑓(𝑡, 𝑥𝑑 , 𝑥𝑎 , 𝑢, 𝑤𝐼 ) 𝑑𝑡 𝑥𝑑 (𝑡𝑘 ) = 𝑥̂𝑑0 𝑐(𝑡, 𝑥𝑑 , 𝑥𝑎 , 𝑢, 𝑤𝐼 ) = 0 ℎ(𝑡, 𝑥𝑑 , 𝑥𝑎 , 𝑢, 𝑤𝐼 ) ≥ 0 𝑦 = 𝑔(𝑡, 𝑥𝑑 , 𝑥𝑎 , 𝑤𝑂 )

(1b) (1c) (1d) (1e) (1f)

Note that Eq. (1c) represents the initial condition of the differential states which is updated at each sampling instant. In the presence of model plant mismatch the simple approach of updating the initial conditions with the value of the states observed, if all the differential states can be measured, can produce a steady-state offset that deteriorates the performance of the controller. Thus, an output feedback method is required to obtain an offset free response

10

. Using the

formulation presented in Eq. (1) the output feedback can be done using an observer such as an extended Kalman filter (EKF) or a moving horizon estimator (MHE) which solution provides the initial conditions of the problem (𝑥̂𝑑0 ) and the value of the disturbances (𝑤𝐼 , 𝑤𝑂 ). There are


10

Page 11 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


other simpler approaches for handling the feedback, when all the differential states are measured, in order to eliminate the offset and they are: the conventional feedback and the steady-state target setting

10

. In the first one, a constant bias, calculated as the difference between the output

predicted and the output measured, is added to the right hand side of Eq. (1f) 4. In the second one, the set-points of the outputs are replaced by a reference value which is obtaining by solving a steady-state optimization problem that minimizes the difference between the reference and the set-point subject to the model that considers the bias 10. If the optimization problem given by Eq. (1) is of large-scale and/or has complex nonlinear terms its solution might take an important amount of time introducing delays and even instabilities in the closed loop control. In addition, the effects of the tuning parameter of the controller on the closed loop performance are not easy to predict and they might not be monotonic 2, 11.

3. NMPC tuning algorithm The optimal tuning problem with respect to the controller parameters can be formulated as follows 11: min

𝑁𝑝 ,𝑁𝑢 ,𝑊,𝜏

𝐽(𝑥𝑑 , 𝑥𝑎 )

(2a)

𝑠. 𝑡. 𝑓̃(𝑥𝑑 , 𝑥𝑎 , 𝑢, 𝑤𝐼 ) = 0 𝑢 = 𝑔̃(𝑥𝑑 , 𝑥𝑎 , 𝑁𝑝 , 𝑁𝑢 , 𝑊, 𝜏) ℎ̃(𝑥𝑎 , 𝑥𝑑 ) ≤ 0

(2b) (2c) (2d)

Where the objective function, Eq (2a), is a performance index for the closed loop behavior of the controller. Eq (2b) is a general set of differential algebraic equations that represent the real plant and it is a function of the states, the control law and the disturbances. This allows to tune the controller for a specific type of anticipated or predefined disturbances. Eq (2d) are inequality constraint for the closed loop behavior (e.g. the maximum overshoot allowed). Note that the


11


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 45

objective function of the problem defined by Eq (2) is an indicator of the whole closed loop performance of the process (e.g. the integral absolute error), which is different from the objective of Eq (1) that represents the OCP objective (e.g. a quadratic function of tracking with respect to a specific set-point). Also, the OCP defined in Eq. 1 is solved recurrently at every sampling interval, while the problem of Eq. 2 is only solved once for the definition of the controller tuning parameters. In most cases it is impossible to define exactly the Eq (2b). Moreover, Eq (2c) is the explicit feedback control law expressed as a function of the differential states, the algebraic states and the controller tuning parameters (e.g. prediction horizon 𝑁𝑝 , control horizon 𝑁𝑢 , objective function weights 𝑊, sampling time 𝜏); and it can only be defined in this way for linear model predictive controllers without inequality constraints. If it were possible to formulate explicitly the optimization problem of Eq (2) its solution would be the tuning parameters of the NMPC controller that minimize the specified closed loop performance index. Therefore, there is a need to formulate the optimization problem of Eq (2) in a different way to find an alternative solution to this problem considering its limitations. The algorithm presented hereafter is one alternative for solving this optimization problem for tuning NMPC controllers, but before introducing it in detail it is important to analyze the tuning parameters of the controller and their general effect on the closed loop performance. 

Sampling time: it is the time interval in which a new control signal is feedback to the system. The sampling interval should be small enough to capture adequately the dynamic of the process, but the smaller it is the optimization problem becomes larger 12. It is not adequate as a tuning parameter because it should be fixed based on the equipment installed and the limitations of sensors and actuators. It also impacts directly those parameters that are


12

Page 13 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


specified in sampling intervals like the prediction and control horizons 5. It is important to consider the Shannon’s sampling theorem when selecting the sampling time because it defines an upper bound for this parameter to avoid loss of information during the sampling of a signal, in other words, it states that the sampling frequency should be greater than two times the maximum frequency of the signal

13

. Note that the frequency of interest here can

be referred to the manipulated variables or to the disturbances, being the first one more important because it indicates how fast the actuators can response. Also, there are some empirical guidelines that suggest a sampling time between 10% and 30% of the period of oscillation for oscillatory systems and of the rise time for non-oscillatory systems 14. 

Prediction horizon: it defines the upper limit of the time horizon for the prediction of the process behavior. If the process model is linear, increasing the prediction horizon has a monotonic effect on the closed loop performance improving its stability. It should be selected as large as possible but finite. As the prediction horizon is increased the size and complexity of the dynamic optimization problem does also increase 1.



Control horizon: this is the number of control actions that are allowed during the OCP and it is lower or equal than the prediction horizon. It affects how aggressive the control action is and defines a trade-off between robustness of the controller for large control horizons and easy of computation of the OCP solution for small control horizons

1, 12

. The control inputs

can also be parameterized in a blocking scheme in which the user specifies points on the control horizon that are not computed and are set equal to the previous control action 4. Before introducing the rest of the NMPC tuning parameters it is important to consider in detail the objective function of Eq. (1). Usually this objective function is composed of different objectives (e.g. quadratic tracking of the set-point, penalization of the rate of change of the


13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 45

manipulated variables) and the compromise among these is expressed in a single objective function constructed as a weighted sum, Eq. (3). Also, Eq. (4) shows the most common objective function used in (N)MPC where 𝑊𝑖 are weights for every individual objective (𝜙𝑖 ). In Eq. (4), the first term is a quadratic tracking of the outputs with respect to a set-point, the first summation is done over the number of outputs 𝑛𝑦 and the second one over the discrete times until the prediction horizon 𝑁𝑝 , and the second term is a penalization of the rate of change of the manipulated variables, the first summation is done over the number of manipulated variables 𝑛𝑢 and the second one over the discrete times until the control horizon 𝑁𝑢 . The second term can also be expressed as the change of the manipulated variables with respect to a reference. Note that this is a discrete time representation 𝑛𝑜

(3)

min 𝑂𝐵𝐽 = ∑ 𝑊𝑖 𝜙𝑖 (𝑡, 𝑦, 𝑢)


𝑖=1 𝑛𝑦

𝑁𝑝

min 𝑂𝐵𝐽 = ∑ 𝑊𝑖 ∑(𝑦𝑖,𝑘 −




𝑖=1

𝑛𝑢 𝑠𝑝 2 𝑦𝑖 )

𝑘=1

𝑁𝑢 −1

+ ∑ 𝑊𝑗 ∑ (Δ𝑢𝑗,𝑘 ) 𝑗=1

2

(4)

𝑘=0

Weights on the outputs: this weights are used in the objective function to scale the variables and to prioritize the different objectives according to their relative importance. Increasing a specific weight, while keeping the rest constant, drives the control effort towards a more tighter control in that variable 1.



Weights on the change or rate of change of the inputs: they suppress the aggressive changes of the manipulated variables making a more robust controller, but they also make the response on the controlled variables slower 5. Setting a small penalty on this part of the objective function produces a more aggressive response and decreases the ability of the controller to handle drastic changes or disturbances.


14

Page 15 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Algorithm for tuning NMPC controllers An algorithm for tuning NMPC controllers is developed in which the final result is the value of all the tuning parameters with the exception of the sampling time that is considered as an input for the algorithm. Figure 2 and Figure 3 present every step of the algorithm and the shadowed blocks are the stages in which the tuning parameter values are obtained. The algorithm consists of three phases. The first one is the definition of the model, the plant and the initialization of some parameters and sets. The second one is the definition of the OCP objective function weights using multiobjective optimization and the computation of key performance indexes (KPI) for each configuration. Note that a KPI can be any type of metric used to quantify the closed-loop performance of the process and some examples are: the integral square error (ISE), the overshoot percentage, the maximum solution time of the OCP and the economic profit. Finally, the third phase is the solution of an optimization problem formulated based on the closed loop performance and the KPIs obtained. Each step of the algorithm is describe next.


15


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 45

PHASE II. Weight definiti construction of KPI set assoc each point in the horizon con

PHASE I. Initialization of model, sets and parameters

k=1

𝐶ℎ𝑜𝑜𝑠𝑒 (𝑁𝑝 , 𝑁𝑢 ) = (𝑁𝑝 , 𝑁𝑢 )𝑘 ; 𝑤ℎ𝑒𝑟𝑒 (

ALGORITHM FOR TUNING NMPC CONTROLLERS

j=1 Define process model, constraints, and real plant for simulation.

𝑂𝐵𝐽 = 𝜙𝑗

Define the state estimation problem and fix an alternative to handle model-plant mismatch

(𝑥𝑑∗ , 𝑥𝑎∗ , 𝑦 ∗ , 𝑢∗ ) = 𝑂𝑝𝑡𝑖𝑚𝑎𝑙 𝑆𝑜𝑙𝑛. 𝑜𝑓 𝑁 𝑤𝑖𝑡ℎ 𝑂𝐵𝐽 = 𝜙𝑗

∗ 𝐶𝑜𝑚𝑝𝑢𝑡𝑒 𝜙𝑖,𝑗 ,𝑘 ∀ 𝑖 ∈ {1, 2, …

Set point tracking

Type of problem

Disturbance rejection More objectives?

Yes. j += 1 Define the desiered set point. If there are set-point changes include them.

Define the disturbances and their behavior over time.

No 𝑈 ∗ 𝐿 ∗ 𝜙𝑖,𝑘 = max(𝜙𝑖,𝑗 ,𝑘 ) , 𝜙𝑖,𝑘 = min(𝜙𝑖,𝑗 ,𝑘 ) 𝑗

𝑗

𝑊𝑖,𝑘 =

Set the sampling time. Consider the Shannon s theorem as an upper limit for this parameter.

1 𝑈 𝐿 ∀𝑖 ∈ {1, 2, 𝜙𝑖,𝑘 − 𝜙𝑖,𝑘 𝑛𝑜

𝑂𝐵𝐽 = ∑ 𝑊𝑖,𝑘 𝜙𝑖

Control objectives set 𝝓 = 𝜙1 , 𝜙2 , … , 𝜙𝑛 𝑜

𝑖=1

(𝑥𝑑∗ , 𝑥𝑎∗ , 𝑦 ∗ , 𝑢∗ ) = 𝑂𝑝𝑡𝑖𝑚𝑎𝑙 𝑆𝑜𝑙𝑛. 𝑜𝑓 𝑁

KPIs set 𝐾𝑃𝐼_𝑆𝑒𝑡 = 𝑂𝐶𝑃 𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛 𝑡𝑖𝑚𝑒, 𝐾𝑃𝐼1 , 𝐾𝑃𝐼2 , … , 𝐾𝑃𝐼𝑛 𝑘

𝑛𝑜

𝑤𝑖𝑡ℎ 𝑂𝐵𝐽 = ∑ 𝑊𝑖,𝑘 𝜙 𝑖=1

Horizon configurations 𝐻_𝑆𝑒𝑡 = (𝑁𝑝 , 𝑁𝑢 )1 , (𝑁𝑝 , 𝑁𝑢 )2 , … , (𝑁𝑝 , 𝑁𝑢 )𝑛

𝐻

𝐶𝑎𝑙𝑐𝑢𝑙𝑎𝑡𝑒 𝑎𝑛𝑑 𝑠𝑎𝑣𝑒 𝐾𝑃𝐼𝑠 𝑓𝑜𝑟

More (Np, Nu) points?

Yes. k += 1

Figure 2. Phase I of the algorithm for tuning NMPC controllers


16

Page 17 of 45

1 2 3 4 5 6 7 8 9 10 ERS11 12 13 al 14 15 an16 ch 17 18 19 urbance 20 rejection 21 22 23Define the disturbances and 24 their behavior 25over time. 26 27 28 29 orem as an 30 31 32 𝑛𝑜 33 34 35 2 , … , 𝐾𝑃𝐼𝑛 𝑘 36 37 , …38 , (𝑁𝑝 , 𝑁𝑢 )𝑛 𝐻 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


PHASE II. Weight definition and construction of KPI set associated with each point in the horizon configuration

PHASE III. Approximation of KPIs and definition of horizons

k=1 𝑓𝑖 (𝑎1 , 𝑎2 , … , 𝑎𝑚 )∗ = 𝑂𝑝𝑡𝑖𝑚𝑎𝑙 𝑆𝑜𝑙𝑛. 𝑜𝑓 𝑅𝑆 − 𝐾𝑃𝐼 (𝑖) 𝐸𝑞 (8) , ∀ 𝑖 ∈ 𝐾𝑃𝐼_𝑆𝑒𝑡

𝐶ℎ𝑜𝑜𝑠𝑒 (𝑁𝑝 , 𝑁𝑢 ) = (𝑁𝑝 , 𝑁𝑢 )𝑘 ; 𝑤ℎ𝑒𝑟𝑒 (𝑁𝑝 , 𝑁𝑢 )𝑘 ∈ 𝐻_𝑆𝑒𝑡

j=1

∗

(𝑁𝑝 , 𝑁𝑢 ) = 𝑂𝑝𝑡𝑖𝑚𝑎𝑙 𝑆𝑜𝑙𝑛. 𝑜𝑓 𝐻𝑃 [𝐸𝑞 (10)] 𝑂𝐵𝐽 = 𝜙𝑗

(𝑥𝑑∗ , 𝑥𝑎∗ , 𝑦 ∗ , 𝑢∗ ) = 𝑂𝑝𝑡𝑖𝑚𝑎𝑙 𝑆𝑜𝑙𝑛. 𝑜𝑓 𝑁𝑀𝑃𝐶 𝐸𝑞 (1) 𝑤𝑖𝑡ℎ 𝑂𝐵𝐽 = 𝜙𝑗

Is (Np, Nu)* in H_set?

∗ 𝐶𝑜𝑚𝑝𝑢𝑡𝑒 𝜙𝑖,𝑗 ,𝑘 ∀ 𝑖 ∈ {1, 2, … , 𝑛𝑜 }

No

Yes

k += 1 More objectives?

Yes. j += 1

∗

𝐻_𝑆𝑒𝑡 ∪ (𝑁𝑝 , 𝑁𝑢 )

END

No 𝑈 ∗ 𝐿 ∗ 𝜙𝑖,𝑘 = max(𝜙𝑖,𝑗 ,𝑘 ) , 𝜙𝑖,𝑘 = min(𝜙𝑖,𝑗 ,𝑘 ) , ∀𝑖 ∈ {1, 2, … , 𝑛𝑜 } 𝑡,𝑗

𝑡,𝑗

𝑊𝑖,𝑘 =

1 𝑈 𝐿 ∀𝑖 ∈ {1, 2, … , 𝑛𝑜 } 𝜙𝑖,𝑘 − 𝜙𝑖,𝑘 𝑛𝑜

𝑂𝐵𝐽 = ∑ 𝑊𝑖,𝑘 𝜙𝑖 𝑖=1

(𝑥𝑑∗ , 𝑥𝑎∗ , 𝑦 ∗ , 𝑢∗ ) = 𝑂𝑝𝑡𝑖𝑚𝑎𝑙 𝑆𝑜𝑙𝑛. 𝑜𝑓 𝑁𝑀𝑃𝐶 𝐸𝑞 (1) 𝑛𝑜

𝑤𝑖𝑡ℎ 𝑂𝐵𝐽 = ∑ 𝑊𝑖,𝑘 𝜙𝑖 𝑖=1

𝐶𝑎𝑙𝑐𝑢𝑙𝑎𝑡𝑒 𝑎𝑛𝑑 𝑠𝑎𝑣𝑒 𝐾𝑃𝐼𝑠 𝑓𝑜𝑟 (𝑁𝑝 , 𝑁𝑢 )𝑘

Yes. k += 1

More (Np, Nu) points?

No

Figure 3. Phase II and III of the algorithm for tuning NMPC controllers In Phase I of the algorithm the model, the real plant and the initialization parameters are defined. The algorithm starts by defining the model of the process. For example, it can be a set of differential equations that describes the dynamic behavior of the system including inequality constraints. These two elements correspond to Eq (1b) – Eq (1e) of the OCP formulation. Within


17


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 45

the tuning algorithm it is important to consider the model plant mismatch because every practical implementation of NMPC exhibits a mismatch and the way this is handled has a significant effect on the tuning parameters of the controller

15

. Therefore, it is necessary to define a

simulation model for the real plant if the tuning cannot be done using the plant itself. Some alternatives for modeling the real plant behavior are: using a more detailed model than the one used for prediction, introducing white noise in the measurements or modifying the parameters of the prediction model. Once this two models, or the plant measurements and the prediction model, are defined the next step is to consider the state estimation problem and an alternative to handle model plant mismatch. In this paper we consider, without loss of generality, that all the states can be measured, thus we do not solve a state estimation problem; but estimators like the Kalman filter, extended Kalman filter or MHE (Moving horizon estimation) can be used in this step of the tuning algorithm

16

. For handling the model plant mismatch and to avoid computational

delays we used the asNMPC strategy coupled with NLP sensitivity which is a fast alternative to update the value of the manipulated variables once the real measurements of the process are available 17. This is described in more detail in Section 4 because it is applied for the case studies presented and the algorithm is not restricted to using it. Also, other alternatives can be used here like a proportional correction of the manipulated variables based on the error in the measured variables 18. The rest of Phase I includes defining the type of control problem for which the controller is being designed (e.g. disturbance rejection or set-point tracking) and initializing some parameters and sets. The parameters that must be initialized are: the sampling time (𝜏), that should be selected based on the system limitations and on the Shannon’s sampling theorem; the set of control objectives (𝝓), which corresponds to each element presented in the weighted sum


18

Page 19 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


objective function, Eq. (3); a set of performance indexes (e.g. integral square error, overshoot percentage, economic profit) (𝐾𝑃𝐼_𝑆𝑒𝑡) and a set of different configurations of prediction horizon and control horizon (𝐻_𝑆𝑒𝑡). Note that the 𝐻_𝑆𝑒𝑡 elements are only an initial configurations of 𝑁𝑝 and 𝑁𝑢 and later in the algorithm more elements will be added to the set. However, the points of this initial configuration must have a good dispersion allowing a good sampling of the whole region defined by 𝑁𝑝 and 𝑁𝑢 in order to fit response surfaces in this space 19

. One alternative to initialize the 𝐻_𝑆𝑒𝑡 is presented in Figure 4, this configuration is based on a

central composed design of experiments for a non-uniform region and allows a good sampling of 19

the space defined by 𝑁𝑝 and 𝑁𝑢

. Note that the sampling space shown in Figure 4 depends on

the value of the parameter 𝑁𝑝𝑚𝑎𝑥 which is an upper bound for the prediction horizon and it is defined by the analyst. Nu

Np = Nu

Np max

(½)Np max

(¼)Np max

2 2

(½)Np max

(¾)Np max

Np max

Np

Figure 4. Initial configuration of H_Set†. In this step of the tuning algorithm additional parameters might be included such as an input blocking scheme or internal model control (IMC) filter constants, but it will increase the number †

For the examples presented later in this paper 𝑁𝑝𝑚𝑎𝑥 takes values of 10 and 50 which means that the sampling is done at values of 𝑁𝑝 equals to (2, 5, 8, 10) and (2, 25, 38, 50), respectively.


19


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 45

of cases that must be evaluated in the Phase II. In other words, the H_Set will be composed not just by the prediction and control horizon, but it will have additional elements corresponding to the additional aspects of the controller. For example, if an input blocking is included in the controller each element of the H_Set will have the form of (𝑁𝑝 , 𝑁𝑢 , 𝐵𝑠), where 𝐵𝑠 is the blocking scheme. Also, adding more items to each element of the H_Set increases the complexity of the Phase III, but it is not a restriction for applying the tuning algorithm proposed. In Phase II of the algorithm for each element of the set 𝐻_𝑆𝑒𝑡, which means for each configuration of prediction horizon and control horizon, the following is done: 

Solve the NMPC problem for each individual control objective (e.g. quadratic tracking of one state variable with respect to its set-point, quadratic penalization of the rate of change of a manipulated variable). At every period of the NMPC, when an OCP is solved, the values ∗ of all the objectives are computed and stored (𝜙𝑖,𝑗,𝑘 ). Also, the notation here indicates that

the objective 𝑖 is computed for every time period when the objective 𝑗 is minimized and the horizon configuration 𝑘 is used. It is important to mention that for a specific set of indexes ∗ (e.g. 𝑖 = 1, 𝑗 = 1, 𝑘 = 1) the variable 𝜙𝑖,𝑗,𝑘 is actually a vector indexed in the discrete time t

because the objective function value is computed at every sampling time. The * symbol indicates the value of the variable/expression at the optimal solution. 

𝑈 Once the NMPC has been solve for every individual objective (𝜙𝑗 ) an upper bound (𝜙𝑖,𝑘 ) 𝐿 and a lower bound (𝜙𝑖,𝑘 ) for each objective are computed from the values stored previously.

Note that the upper and lower bounds are computed using the maximum and minimum ∗ operators over the 𝑗 index and the discrete time (𝑡) on the variable 𝜙𝑖,𝑗,𝑘 . Then, using these

bounds the range of all objectives is computed. 

The OCP objective function weights are calculated based on the range of each objective.


20

Page 21 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




Another NMPC problem is solved using a composed objective function which is defined as the weighted sum of all the objectives and for this final result the KPIs are computed and saved. The las three steps of Phase II of the proposed algorithm are better explained in the following

𝑈 𝐿 paragraphs and in the Appendix A which shows in great detail how to compute 𝜙𝑖,𝑘 and 𝜙𝑖,𝑘 for

the case study 1 (Section 5.1). At the end of Phase II a set of weights and a set of KPIs must have been defined for each configuration of prediction and control horizon. The idea behind this way of calculating the objective function weights is a multiobjecive optimization technique called utopia tracking approach which minimizes the distance between the Pareto front and the utopia point 20. Figure 5 illustrates the utopia point, the nadir point, the Pareto font and the distance that is minimized using this approach for a case consisting only of two objectives. Zavala and Flores use this utopia tracking approach to solve online the multiple conflicting objectives in a NMPC controller. They solve in parallel as many optimization problems as objectives are defined to find the utopia point and then the OCP that is solved to define the feedback value of the manipulated variables has the following objective function 21: 𝑛𝑜

min ∑ |


𝑖=1

𝜙𝑖 − 𝜙𝑖𝐿 | 𝜙𝑖𝑈 − 𝜙𝑖𝐿

(5)

Note that in Eq. (5) is expressed in terms of the elements of the following vectors: 𝝓 is a vector of objectives, 𝝓𝑳 is a vector of lower bounds for the objectives (utopia point) and 𝝓𝑼 is a vector of upper bounds for the objectives (nadir point). Also, each distance from the objective to the utopia point is standardized dividing it by the range of the objective because the order of


21


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 45

magnitude of the objectives might differ significantly and produces an ill-conditioned optimization problem. 𝜙2𝑈

𝜙1𝑈

Nadir point

Pareto front

𝜙1 Compromise solution

Utopia point

𝜙1𝐿

𝜙2𝐿

𝜙2

Figure 5. Illustration of the utopia tracking approach for two objectives. It is possible to avoid the online solution of more than one optimization problem within every time interval if the global coordinates of the utopia and the nadir point are defined previously. In order to compute the utopia and the nadir point it is necessary to solve a NMPC problem for each objective and, once this is done, it is possible to use Eq (5) as the objective function of the OCP. However, considering that the utopia point is an unreachable point because it minimizes all the objectives simultaneously

22

, the absolute value operator of Eq (5) can be omitted and the

objective function can be rearranged as follows: 𝑛

𝑛

𝑛

𝑛

𝜙𝑖 − 𝜙𝑖𝐿 𝜙𝑖 − 𝜙𝑖𝐿 𝜙𝑖 𝜙𝑖𝐿 min ∑ | 𝑈 | = min ∑ = min ∑ − ∑ 𝑥𝑑 ,𝑥𝑎 ,𝑢 𝑥𝑑 ,𝑥𝑎 ,𝑢 𝜙𝑖 − 𝜙𝑖𝐿 𝜙𝑖𝑈 − 𝜙𝑖𝐿 𝑥𝑑 ,𝑥𝑎,𝑢 𝜙𝑖𝑈 − 𝜙𝑖𝐿 𝜙𝑖𝑈 − 𝜙𝑖𝐿 𝑖=1

𝑖=1

𝑖=1

(6)

𝑖=1

Finally, taking into account that the argument of minimizing a function plus a constant is the same as minimizing only the function (it is also valid if the problem is constrained)

23

the OCP


22

Page 23 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


objective function can be expressed as a weighted sum of the individual objectives where each weight is the inverse of the objective range, Eq (7). (𝑥𝑑∗ , 𝑥𝑎∗ , 𝑢∗ ) =

𝑛

𝑛

𝑛

𝑛

𝑖=1

𝑖=1

𝑖=1

𝑖=1

𝜙𝑖 𝜙𝑖𝐿 𝜙𝑖 arg min arg min arg min ∑ 𝑈 − ∑ = ∑ 𝑈 +𝐶 = ∑ 𝑊𝑖 𝜙𝑖 𝑈 𝐿 𝐿 (𝑥𝑑 , 𝑥𝑎 , 𝑢) (𝑥𝑑 , 𝑥𝑎 , 𝑢) (𝑥𝑑 , 𝑥𝑎 , 𝑢) 𝜙𝑖 − 𝜙𝑖 𝜙𝑖 − 𝜙𝑖 𝜙𝑖 − 𝜙𝑖𝐿

(7)

When the utopia tracking approach is used for defining the OCP weights all the objectives receive the same relative importance 24. The Phase III of the algorithm is composed of two steps. The first one is the computation of a response surface that is a function of the prediction horizon (𝑁𝑝 ) and the control horizon (𝑁𝑢 ) for every KPI in the 𝐾𝑃𝐼_𝑆𝑒𝑡. Note in Figure 4 that only specific measurements of the KPI are available after Phase II of the proposed tuning algorithm and there is not an explicit function of the KPI in the whole space defined by 𝑁𝑝 and 𝑁𝑢 . Because of this it is necessary to use these ‘measured’ points, which are computed after each iteration of Phase II, to approximate a function or response surface that can predict accurately the behavior of the KPIs. To compute the response surface it is necessary to solve the problem RS-KPI (Response Surface for KPI), stated in Eq (8). This problem minimizes the deviation between the ‘measured’ points (obtained in Phase II) and a specific function; in other words, it is a linear or nonlinear fit. In this paper we use a quadratic approximation as it is shown in Eq (9) to describe how each KPI depends on 𝑁𝑝 and 𝑁𝑢 , but other types of response surfaces can be used for this approximation. min

𝑎1 ,𝑎2 ,…,𝑎𝑚

̂𝑖 − 𝐾𝑖 )2 ∑ (𝐾

(8)

𝑖∈𝐻_𝑆𝑒𝑡

𝑠. 𝑡. 𝐾𝑖 = 𝑓(𝑎1 , 𝑎2 , … , 𝑎𝑚 ), ∀𝑖 ∈ 𝐻_𝑆𝑒𝑡 𝑓(𝑎1 , 𝑎2 , … , 𝑎𝑚 ) = 𝑎20 𝑁𝑝2 + 𝑎10 𝑁𝑝 + 𝑎02 𝑁𝑢2 + 𝑎01 𝑁𝑢 + 𝑎11 𝑁𝑝 𝑁𝑢 + 𝑎00

(9)

The second step of the Phase III is to solve the horizons problem (HP) which is similar to the problem stated by Eq (2), but the only decision variables are 𝑁𝑝 and 𝑁𝑢 because the objective


23


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 45

function weights have been defined previously. The HP minimizes a KPI of the closed loop behavior of the system subject to the nature of the control and prediction horizon, to a constraint in the solution time of the OCP, and to other inequality constraints based on the performance indexes (e.g. maximum allowable overshoot percentage, minimum rise time). Eq (10) presents a general formulation of the HP, where the function 𝐹 is a KPI function obtained from solving the RS-KPI problem, Eq. 8, for a specific performance index. G is a vector of KPI functions and each one corresponds to the solution of a RS-KPI problem. min 𝐹(𝑁𝑝 , 𝑁𝑢 )

𝑁𝑝 ,𝑁𝑢

𝑠. 𝑡. 𝑡𝑠𝑜𝑙 ≤ 𝜏 𝐺(𝑁𝑝 , 𝑁𝑢 ) ≤ 0 𝑁𝑢 ≤ 𝑁𝑝 ≤ 𝑁𝑝𝑚𝑎𝑥 (𝑁𝑢 , 𝑁𝑝 ) ∈ ℤ+

(10)

The final step of the algorithm is to check if the solution of the HP is already in the 𝐻_𝑆𝑒𝑡, if so the algorithm finish and the arguments of the solution of the HP are the prediction and control horizon for the NMPC controller. Note that this pair of points have a vector of objective function weights associated which were found in Phase II. If the solution of HP is not in the 𝐻_𝑆𝑒𝑡 this new point is added to the set and the Phase II is repeated for this point to define the corresponding objective function weights.

4. NLP sensitivity and asNMPC controller The ideal NMPC assumes that the OCP can be solved instantaneously, that all the states are measured and that there is no model plant mismatch. However, if all these assumptions are considered valid during the tuning of the controller the final result might not be robust and the controller could fail even under small disturbances. Thus, it is necessary to consider an additional


24

Page 25 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


level of complexity. For doing so, in this paper we consider the advanced step NMPC (asNMPC) controller that avoids computational delays and it is capable of handling model plant mismatch. The asNMPC solves the OCP in advanced, which means that at the current time 𝑡𝑘 and with the measurements of the current states and manipulated variables, the prediction model is used to predict the states at the next time period 𝑡𝑘+1 and the OCP for this last period is solved between 𝑡𝑘 and 𝑡𝑘+1 . Using this strategy the delays generated due to the computational time are avoided, but there can be still a mismatch between the states predicted and the measured values at 𝑡𝑘+1 , thus a rapid correction using NLP sensitivity has to be done. If the parameters (𝑝) of an optimization problem, which in the case of the OCP correspond to the initial conditions and disturbances, are associated with artificial constraints of the form 𝑥𝑝 − 𝑝 = 0 where 𝑥𝑝 is a variable, the NLP sensitivity calculations can be reduced to solve the linear system given by Eq. (11)

25

. Note that this is a linear equation in which 𝐻𝐿 is the Hessian matrix of the Lagrange

function (𝐿), 𝐴 is the Jacobian matrix of the equality constraints, 𝑉 is a diagonal matrix which contains the dual variables (𝑣), 𝐿 is the Lagrange function, 𝑐 is a vector of equality constraints, 𝑍 is a diagonal matrix which contains all the variables of the problem and 𝜆 are the Lagrange multipliers for the equality constraints. Eq. (11) is obtained after applying the chain rule over the first order optimality conditions of a constrained NLP problem when it is solved using an interior point algorithm, such as the one implemented in the solver IPOPT. These derivatives are taken with respect to the parameters in order to obtain the sensitivity of the variables (𝑧, 𝑥𝑝 ). Also, the left hand side matrix corresponds to the KKT matrix at the optimal solution which is available after solving the optimization problem. For a more detailed discussion about asNMPC controller, its properties and NLP sensitivity calculations the reader is referred to the following references 25, 26, 27

.


25


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

∇𝑧𝑥𝑝 𝐿(𝑧, 𝑥𝑝 , 𝜆, 𝑣)

𝐴

−𝐼

∇𝑥𝑝 𝑧 𝐿(𝑧, 𝑥𝑝 , 𝜆, 𝑣)

∇𝑥𝑝 𝑥𝑝 𝐿(𝑧, 𝑥𝑝 , 𝜆, 𝑣)

∇𝑥𝑝 𝑐(𝑧, 𝑥𝑝 )

0

0 0 0

0 𝑍 0

[

𝐻𝐿 𝐴𝑇 𝑉 0

∇𝑥𝑝 𝑐(𝑧, 𝑥𝑝 ) 0 𝐼

𝑇

Page 26 of 45

0

Δ𝑧 0 𝐼 Δ𝑥𝑝 0 0 = 0 Δ𝜆 0 Δ𝑣 0 [ Δ𝜆̅ ] [Δ𝑝] 0]

(11)

Along this manuscript we focus on NMPC and thus the OCP, after discretization of the differential equations, is defined as a large-scale NLP. The solver IPOPT is chosen for solving this NLP problem because it has been proven effective for the type of formulations that arise in chemical processes

23, 28

and it has been apply in NMPC controllers with good results

29

. The

IPOPT solver implements a deterministic optimization algorithm based on the Newton’s method for solving the KKT conditions of the problem. Moreover, it uses a primal-dual interior point algorithm, in which the inequality constraints are handled in the objective function introducing a logarithm barrier and a barrier parameter, and a filter line search method for solving the first order optimality conditions. sIPOPT is a toolbox that can be found as an extension of IPOPT solver and that performs the NLP sensitivity calculations automatically, Eq. (11). The only requirement for the user is to define the artificial constraints within the problem formulation and a series of suffixes used to communicate the problem with the solver indicating which of the variables are artificial. For further information about the solvers IPOPT and its extension sIPOPT the reader is referred to the references 25, 30.

5. Examples This section presents three examples that show the application of the tuning algorithm described in Section 3. All the optimization problems were formulated in Pyomo

31

, an

optimization software imbedded in Python, and solved with the nonlinear solver IPOPT, the NLP


26

Page 27 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


sensitivity calculation was performed with the extension of this solver sIPOPT. The whole tuning algorithm was programed in Python. 5.1.

Case 1: CSTR with two control objectives

This is an example adapted from

32

in which a CSTR is used for the production of propylene

glycol and it is desired to control the reactor temperature during the start-up operation and in the presence of disturbances. The hydrolysis reaction of propylene oxide to produce propylene glycol is carried out with excess of water, in the presence of methanol and using a strong acid catalyst. The reaction is first-order in propylene oxide concentration. The dynamic model and operational constraints for this example are given by Eq (12) where A: Propylene oxide, B: Water, C: Propylene glycol, M: Methanol. This model is based on molar balances Eq (12a) and energy balances Eq (12b) assuming that the hold-up of cooling fluid in the jacket is negligible; also the parameters of the model are presented in Table 1. Here 𝐶𝑖 is the molar concentration of specie 𝑖, 𝑟 is the reaction rate, 𝑉 is the reactor volume, 𝜗𝑖 is the stoichiometric coefficient of species 𝑖, 𝑣0 is the volumetric feed flow rate, 𝑇 is the reactor temperature, 𝐶𝑝 is the heat capacity, 𝐹 is the molar flow rate, Δ𝐻𝑟𝑥𝑛 is the reaction enthalpy, 𝑇𝑤 is the inlet temperature of the coolant fluid, 𝑈 is the global heat transfer coefficient and 𝐴 is the heat transfer area. The objective of this control problem is to maintain an operational temperature of 63°C manipulating the coolant fluid feed flow rate (𝑚̇𝑤 ). Note that the reactor do not start from its steady state operation point, thus it is a start-up operation and after two hours a disturbance is introduced; more specifically the inlet temperature drops linearly from 24°C to 21°C in an one hour interval.


27


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

𝑑𝐶𝑖 𝑣0 = 𝜗𝑖 𝑟 + (𝐶𝑖0 − 𝐶𝑖 ) , ∀𝑖 ∈ {𝐴, 𝐵, 𝐶, 𝑀} 𝑑𝑡 𝑉 𝑑𝑇 𝑄 − ∑𝑖∈{𝐴,𝐵,𝐶,𝑀} 𝐹𝑖0 𝐶𝑝,𝑖 (𝑇 − 𝑇0 ) − (Δ𝐻𝑟𝑥𝑛 )𝑟𝑉 = ∑𝑖∈{𝐴,𝐵,𝐶,𝑀} 𝐶𝑖 𝑉𝐶𝑝,𝑖 𝑑𝑡 𝐸 𝑟 = 𝑘0 exp − 𝐶 𝑅𝑇 𝐴 𝑈𝐴 𝑄 = 𝑚̇𝑤 𝐶𝑝,𝑤 (𝑇𝑤,𝑖𝑛 − 𝑇) [1 − exp (− )] 𝑚̇𝑤 𝐶𝑝,𝑤 𝑇 ≤ 82°𝐶

Page 28 of 45

(12a) (12b) (12c) (12d) (12e)

Table 1. Model parameters for case 1 and 2

𝑉

Parameter

Value 1.9 𝑚3

𝑣0

12.5 𝑚3 /ℎ

𝐹𝑖0 {𝐴, 𝐵, 𝐶, 𝑀} 𝐶𝑝,𝑖 {𝐴, 𝐵, 𝐶, 𝑀} 𝑇0 𝐶𝑖 (𝑡 = 0){𝐴, 𝐵, 𝐶, 𝑀}

𝑘𝑚𝑜𝑙 ℎ 𝑘𝐽 {146.5, 75.2, 192.6, 81.72} 𝑘𝑚𝑜𝑙. 𝐾 24°𝐶 𝑘𝑚𝑜𝑙 {0, 4.64, 0, 0} 𝑚3 {36.3, 453.6, 0, 45.36}

Parameter 𝑘0 𝐸 𝐶𝑝,𝑤 𝑈𝐴 𝑇𝑤,𝑖𝑛 Δ𝐻𝑟𝑥𝑛

Value 16.96𝑥1012 ℎ−1 𝑘𝐽 75362.4 𝑘𝑚𝑜𝑙 𝑘𝐽 75.3 𝑘𝑚𝑜𝑙. 𝐾 𝑘𝐽 30386 ℎ. 𝐾 15.5°𝐶 𝑘𝐽 −83736 𝑘𝑚𝑜𝑙

The real plant is considered to be the same model described by Eq. 12 but decreasing the frequency factor of the reaction rate expression by 2%. Additional model plant mismatch is introduced adding white noise of amplitude 1% to the measurements of the concentration and temperature at each sampling period. According to the Phase I of the algorithm proposed in Section 3 after defining the type of control problem, the real plant model for simulation and the prediction model it is necessary to define the sampling time, the control objectives and the KPIs that evaluate the closed-loop performance of the controller. For this example a sampling time of 0.01 h is selected which is small enough considering that the time constant of this process for changes around its final steady state is approximately 0.09 h. Two control objectives are defined for this problem: a quadratic tracking of the temperature set-point within the prediction horizon, Eq. (13), and a


28

Page 29 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


penalization of the drastic changes of the coolant feed flow rate within the control horizon, Eq. (14). Finally, three KPIs are selected to evaluate the closed-loop behavior of the process: the integral square error (ISE) of the temperature tracking, the maximum solution time of the OCP within the NMPC scheme and the average concentration of propylene glycol (C) during the time window of the operation. 𝜙1 = ∑

𝑁𝑝

(13)

(𝑇𝑖 − 𝑇 𝑠𝑝 )2

𝑖=0 𝑁𝑢 −1

(14)

𝜙2 = ∑ (Δ𝑚̇𝑤 )2 𝑖=0

Using these definitions the Phase I of the algorithm is over and one can proceed to execute Phase II and III. Nevertheless, before doing so it is necessary to define the horizons problem (HP) which in this case is the minimization of the ISE of the temperature tracking, Eq. (15a), ̅̅̅𝐶 ), Eq. (15b), to a subject to a lower limit in the average concentration of propylene glycol (𝐶 upper limit in the time allowed for solving (𝑡𝑠𝑜𝑙 ) the OCP, Eq (15c), and to the nature of the prediction and control horizons, Eq (15d). Note that every element of the optimization problem defined in Eq (15) are KPIs of the NMPC closed-loop performance and a response surface is defined for each one after solving the RS-KPI problem, in this case we use a quadratic approximation for each function. Examples of these response surface functions are presented in Figure 6 which are constructed for the final iteration of the algorithm and from top to bottom is the OCP maximum solution time, the average concentration of propylene glycol (C) and the ISE of the temperature tracking. In addition, the quadratic fitting of the response surfaces is good enough for these cases producing a 𝑅 2 statistic greater than 0.9 and, furthermore, this functions capture the behavior of the KPIs as a function of the horizons which is not linear nor monotonic. 𝑡=𝑡𝑓

min ∫

2

(15a)

(𝑇 − 𝑇𝑠𝑝 ) 𝑑𝑡

𝑡=0


29


𝑘𝑚𝑜𝑙 𝑚3 𝑡𝑠𝑜𝑙 ≤ 𝜏 = 0.01 ℎ 𝑁𝑢 ≤ 𝑁𝑝 ≤ 𝑁𝑝𝑚𝑎𝑥 = 50

(15b)

𝑠. 𝑡. ̅̅̅ 𝐶𝐶 ≥ 2.35

OCP solution time [s]

(15c) (15d)

6 4 2 0 50 40 30 20 10

Average C concentration [kmol/m3]

Nu

0

0

10

20

30

40

50

Np

2.5 2.4 2.3 2.2 2.1 50 40 30 20 10 Nu

x 10

ISE temperature tracking

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 45

0

0

10

20

30

40

50

Np

4

10

5

0 50 40 30 20 10 Nu

0

0

10

20

30

40

50

Np

Figure 6. Response surface for the KPIs in the case 1


30

Page 31 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


After taking all these definition into account and executing Phase II and III of the algorithm the final tuning parameters are: 𝑁𝑝 = 𝑁𝑢 = 11, 𝑊1 = 2.64𝑥10−5, and 𝑊2 = 4.32𝑥10−4. For a step by step computation of the objective function weights the reader is referred to the Appendix A. First of all, note that the final configuration of prediction horizon and control horizon is not included in the initialization of 𝐻_𝑆𝑒𝑡, this is a result of solving the HP problem every time the Phase III of the algorithm is executed. In addition, with every iteration of Phase III the response surfaces of the KPIs are improved by adding a new point to the sampling space, and this new point is closer to the optimal solution of the HP. Second of all, the weights for the composed objective function are the result of Phase II that executes the utopia tracking multiobjective approach for the configuration 𝑁𝑝 = 𝑁𝑢 = 11 and scale all the objectives to have the same level of importance. 5.2.

Case 2: CSTR with three control objectives

This example uses the same model and considerations of the case 1 with the only exception that an additional objective for the OCP is introduced and it is the maximization of the propylene glycol concentration at the final time of the prediction, Eq. (16). As the utopia tracking approach considers all the objectives in the minimization sense it is necessary to minimize the negative of the concentration in order to maximize it. Introducing this additional objective increases the complexity of the tuning algorithm as an additional weight has to be determined and, during Phase II, an extra NMPC problem has to be solved. However, this shows that the proposed tuning algorithm is capable of handling as many objectives as the analyst decide to include even if it is a nonstandard objective such as Eq. (16) that represents a final cost. (16)

𝜙3 = −𝐶𝑐 (𝑡 = 𝜏𝑁𝑝 )


31


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 45

The results of executing the tuning algorithm with this change are clearly different from those of the case 1 and are the following: 𝑁𝑝 = 19, 𝑁𝑢 = 14, 𝑊1 = 1.63𝑥10−5 , 𝑊2 = 4.42𝑥10−7 and 𝑊3 = 0.376. Adding one extra objective to the problem formulation changes the prediction horizon, the control horizon and the weights of the common objectives since the priorities in the objective function and the velocity with which the controller response to disturbances change. Comparing case 1 and case 2 the value of 𝑊1 does not change significantly but the value of 𝑊2 is reduced allowing more drastic and fast changes of the cooling fluid feed flow rate to achieve a tighter temperature control and thus maximize the concentration of the desired product at the end of the prediction horizon, in fact the ISE of the temperature tracking is improved by a 22.56% in case 2. Moreover, the horizon lengths for case 2 are larger than for case 1 which introduces more control actions during the prediction horizon increasing the robustness of the controller. The profiles of the final tuning configuration for the asNMPC controllers are shown in Figure 7 and it also compares the results of the case 1 and case 2. In addition, Figure 7 shows a tighter temperature control for case 2 than for case 1 that is achieved with faster and more drastic changes of the manipulated variable. Using the tuning parameters obtained from the algorithm proposed in this paper it is possible to maintain the reactor temperature around its set-point even in the presence of noise, disturbances and model plant mismatch without drastic changes in the cooling fluid flow rate.


32

Page 33 of 45

90

a) Case 1 Case 2 sp

Temperature [°C]

80 70 60 50 40 30 20 0

1

2

3

4

5

Time [h] 700

Cooling fluid feed flow rate [kmol/h]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


b) Case 1 Case 2

600 500 400 300 200 100 0 0

1

2

3

4

5

Time [h]

Figure 7. Profiles for the controlled variables (a) and the manipulated variable (b) for case 1 and case 2. 5.3.

Case 3: Extractive distillation column

This final example is a larger and more complex problem than the previous two which serves to show that the NMPC tuning algorithm proposed is capable of handling highly non-linear models and active inequality constraints that produces complex interactions of the tuning parameters of the controller. This case is the design of a NMPC controller with economic orientation for an extractive distillation column used in the separation of the azeotropic mixture of water and ethanol using glycerol as an entrainer. More details about this case and a detail description of the dynamic model used for prediction can be found in a previous work of the authors

33

. In summary, the


33


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 45

model is composed by a set of differential algebraic equations that represent the mass balances, energy balances, phase equilibrium equations and hydraulic relations for each stage of the column. Also, the system is subject to a minimum purity constraint of the distillate, which is a molar fraction of ethanol greater or equal than 99.5%. The Phase I of the algorithm uses the same initialization for the 𝐻_𝑆𝑒𝑡 as the previous two cases and the controller is designed for disturbance rejection using as benchmark disturbance a sinusoidal disturbance in the ethanol feed composition of amplitude 5% and period of 120 min. Additional model plant mismatch is introduced by adding a random noise of amplitude 5% in the liquid hold-up values at each sampling period. The sampling time is chosen as 2.5 min which is close to the time constant of the process for changes around its steady state (≈2 min), but lower than the rise time (≈6 min). Four objectives are considered for the OCP formulation Eq. (17) – (20) among them is one economic objective of profit maximization (it includes the profit from distillate 𝐷 and the cost from the reboiler duty 𝑄𝑅 ), one tracking objective of ethanol mole fraction in distillate (𝑥𝐷𝐸𝑡𝑂𝐻 ), and two penalizations of the changes of the manipulated variables, the reflux ratio (𝑅) and the reboiler duty (𝑄𝑅 ), with respect to a reference value. Finally, three KPIs are selected to evaluate the closed-loop behavior of the process: the total economic profit, the maximum solution time of the OCP within the NMPC scheme and the average ethanol mole fraction in distillate. 𝜙1 = − ∑ 𝜙2 = ∑ 𝜙3 = ∑ 𝜙4 = ∑

𝑁𝑝

(17)

𝐶𝐷 𝐷𝑖 − 𝐶𝑄 𝑄𝑅,𝑖

𝑖=0

𝑁𝑝

𝐸𝑡𝑂𝐻 (𝑥𝐷,𝑖 − 𝑥𝐷𝑠𝑝 )

𝑖=0 𝑁𝑢

2

(𝑄𝑅,𝑖 − 𝑄𝑅,𝑟𝑒𝑓 )

𝑖=0 𝑁𝑢

(𝑅𝑖 − 𝑅𝑟𝑒𝑓 )

(18)

2

(19)

2

(20)

𝑖=0


34

Page 35 of 45

For executing Phase II and III of the algorithm it is necessary to formulate the horizons problem (HP) which is the maximization of the economic profit subject to the ethanol purity constraint in distillate and to the nature of the prediction and control horizon as it is shown in Eq. (21). Every element in the HP is a KPI of the closed loop behavior and thus it is approximated with a response surface solving the RS-KPI problem. Examples of these response surfaces are presented in Figure 8 which are constructed for the final iteration of the algorithm and from top to bottom is the OCP maximum solution time, the average ethanol percentage in distillate and the economic profit. 𝑡=𝑡𝑓

min − ∫

(21a)

(𝐶𝐷 𝐷 − 𝐶𝑄 𝑄𝑅 )𝑑𝑡

𝑡=0 𝑠. 𝑡. ̅̅̅̅̅̅̅̅ 𝑥𝐷𝐸𝑡𝑂𝐻 ≥

(21b) (21c) (21d)

99.5% 𝑡𝑠𝑜𝑙 ≤ 𝜏 = 2.5 𝑚𝑖𝑛 𝑁𝑢 ≤ 𝑁𝑝 ≤ 𝑁𝑝𝑚𝑎𝑥 = 10

OCP solution time [min]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4 3 2 1

0 10 8 6 4 2 Nu

0

0

2

4

6

8

10

Np


35

Page 36 of 45

99.65 99.6 99.55 99.5 10 8 6 4 2 Nu

0

0

2

4

6

8

10

Np

4800

Economic profit ($)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Average ethanol in distillate [%]


4600 4400 4200 10 8 6 4 2 Nu

0

0

2

4

6

8

10

Np

Figure 8. Response surface for the KPIs in the case 3 After taking all these definitions into account and executing Phase II and III of the algorithm the final tuning parameters of the asNMPC controller are: 𝑁𝑝 = 8, 𝑁𝑢 = 5, 𝑊1 = 0.0072, 𝑊2 = 5.03, 𝑊3 = 0.001 and 𝑊4 = 327.91. Note that using this tuning algorithm the order of magnitude of the objective weights can vary significantly due to the differences in scales of every objective. For example, the economic objective, Eq (17), is in the order of 103 , the mole fraction is expressed as a percentage so its order is of 102 , the reboiler duty is of the order of 102 and the reflux ratio is of the order of 10−1. However, this differences in the weights allows a better prioritization of the conflicting objectives presented in the composed objective function of the OCP.


36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 9 shows a comparison of the closed loop profiles using the tuning parameters obtained from the proposed algorithm and a set of arbitrary parameters defined by trial and error, as those defined in a previous work

33

. Using a trial and error approach is difficult to define the

penalization weights for the manipulated variables in order to use all of them to control a specific variable and in this case Figure 9a and Figure 9b show how poorly tuned penalization weights drive all the control effort towards only one manipulated variable, the reflux ratio, even when two variables are available. The proposed tuning algorithm avoids this kind of behavior and defines the composed objective function weights such as all the parameters receive the same level of importance and all the manipulated variables can be used to control the process. Also, the tuning algorithm for this case aims to maximize the economic profit of the process while satisfying the constraints and an increase of the 1.3% is obtained. Note that this improvement in economic profit is obtained even when the prediction and control horizons are smaller than those defined in the trial and error approach. 102

a) Trial and error Tuning algorithm

100

Reboiler duty [MJ/min]

Page 37 of 45

98 96 94 92 90 88 86 0

50

100

150 200 Time[min]

250

300


37


1

b) Trial and error Tuning algorithm

0.9

Reflux ratio [-]

0.8 0.7 0.6 0.5 0.4 0

Ethanol percentage in distillate [%]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 45

100

50

100

150 200 Time[min]

250

300

c) Trial and error Tuning algorithm

99.8 99.6 99.4 99.2 99 0

50

100

150 200 Time[min]

250

300

Figure 9. Profiles for the case 3 and comparison of the asNMPC controller performance using the tuning algorithm versus a trial and error alternative. a) Reboiler duty, b) Reflux ratio, c) Ethanol mole fraction in distillate.

6. Conclusions and perspectives This paper presents a general tuning algorithm for NMPC controllers that consists of three main phases. The first one is the definition of the model and some parameters, the second one is the calculation of the OCP objective function weights using a multiobjective optimization approach based on the utopia tracking concept, and the third one is the definition of the horizons solving an optimization problem that considers the KPIs of the closed loop performance of the controller. Within the tuning algorithm the model plant mismatch is considered and it is handled using the asNMPC scheme coupled with NLP sensitivity, but other strategies to handle the


38

Page 39 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


mismatch problem can be implemented according to the analyst. This tuning algorithm is a versatile approach that can be used for defining systematically the parameters of NMPC controllers regardless the prediction model, the type of constraints, the control objectives or the way the estimation problem and the mismatch problem are handled. The application and effectiveness of the proposed tuning algorithm are illustrated with three examples, two of them are the control of a CSTR reactor during a start-up operation and the third one is the control of an extractive distillation column under the presence of disturbances. These three cases show how to apply the tuning algorithm, the considerations that must be taken into account while applying it and how it can handle nonlinearities, inequality constraints and economic objectives. For the three cases the NMPC tuning parameters defined by the algorithm allow a good closed loop performance of the systems, even in the presence of model plant mismatch cause by random noise or inaccurate predictions. Also, when the OCP has economic objectives or final cost objectives the proposed tuning algorithm define efficiently the weight for these objectives. Even though the examples presented in this paper are focus on the application of the NMPC tuning algorithm to chemical processes, this methodology can be extended to others applications. Some extension can be further incorporated in the algorithm to improve its results such as: the incorporation of the state estimation problem, using a better representation for the response surface of the KPIs, and other multiobjective optimization technique to define the OCP objective function weights.

Appendix A. Detail results and calculations for Case 1 This appendix aims to explain in a more detailed level the Phase II of the proposed algorithm using the case study 1 as an example.


39


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 45

When a horizon configuration is chosen (𝑁𝑝 , 𝑁𝑢 )𝑘 , for the results presented in this appendix is (11, 11), the Phase II deals with the problem of defining the objective function weights using the utopia tracking approach. First of all, the NMPC problem is solved as many times as objectives are with 𝜙𝑗 as objective function and the results for the case study 1, until this step, are presented in Table A 1. Second, the lower and upper bounds of each objective is calculated and these values are used to computes the weights. Note that the values presented at the bottom of the table ∗ are the minimum and maximum values of 𝜙𝑖,𝑗,𝑘 with respect to the discrete time (𝑡), but they do

not correspond the lower and upper bounds of the objectives; an additional calculation is required. Table A 1. Results for solving the NMPC minimizing 𝜙𝑗 for the case study 1 Sampling period

Time [h]

1

Minimizing 𝝓𝟏

Minimizing 𝝓𝟐

∗ 𝜙1,1,𝑘

∗ 𝜙2,1,𝑘

∗ 𝜙1,2,𝑘

∗ 𝜙2,2,𝑘

0

3.79.E+04

0

3.79.E+04

1.21.E-08

2

0.01

3.73.E+04

0

3.73.E+04

8.85.E-09

3

0.02

3.64.E+04

0

3.64.E+04

6.79.E-09

4

0.03

3.51.E+04

0

3.51.E+04

5.47.E-09

5

0.04

3.37.E+04

0

3.37.E+04

4.68.E-09

6

0.05

3.32.E+04

0

3.32.E+04

6.50.E-09

7

0.06

3.23.E+04

0

3.23.E+04

2.61.E-09

8

0.07

3.13.E+04

0

3.13.E+04

1.67.E-09

9

0.08

3.07.E+04

0

3.07.E+04

4.24.E-09

10

0.09

2.97.E+04

0

2.97.E+04

3.38.E-09

⋮

⋮

⋮

⋮

⋮

⋮

498

4.97

1.84.E+00 1.43.E+03 3.28.E+03

2.82.E-21


40

Page 41 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


499

4.98

6.04.E-01

1.01.E+03 3.42.E+03

2.83.E-21

500

4.99

5.16.E-01

4.81.E+02 3.66.E+03

2.89.E-21

501

5.00

6.90.E-04

1.70.E+02 3.30.E+03

2.78.E-21

Min

4.89.E-04

Max

3.79.E+04 2.31.E+03 3.79.E+04 1.98.E+02

0

3.08.E+02

0

𝑈 𝐿 In order to calculate the upper (𝜙𝑖,𝑘 ) and lower (𝜙𝑖,𝑘 ) bounds of the objectives it is necessary ∗ to compute the maximum and minimum of 𝜙𝑖,𝑗,𝑘 with respect to the index 𝑗 and to the discrete

time. In Table A 1 only the minimum and maximum with respect to the discrete time are presented and it is necessary to compute the minimum and maximum with respect to the index 𝑗 for every objective 𝑖. Table A 2 summaries these results for the case study 1. Finally, the weights 𝑈 𝐿 calculated according to Eq. (7) and the values obtained for 𝜙𝑖,𝑘 and for 𝜙𝑖,𝑘 are: 𝑊1 =

2.64𝑥10−5 and 𝑊2 = 4.32𝑥10−4 . Table A 2. Calculation of 𝝓𝑳 and 𝝓𝑼 for case study 1

𝒎𝒊𝒏 𝒕 𝒎𝒂𝒙 𝒕

𝝓∗𝟏,𝒋,𝒌 𝝓∗𝟐,𝒋,𝒌 𝝓∗𝟏,𝒋,𝒌 𝝓∗𝟐,𝒋,𝒌

Minimizing 𝝓𝟏 (𝒋 = 𝟏) 4.89.E-04 0 3.79.E+04 2.31E+03

Minimizing 𝝓𝟐 (𝒋 = 𝟐) 3.08.E+02 0 3.79.E+04 1.98E+02

𝝓𝑳𝒊,𝒌 4.89.E-04 0 X X

𝝓𝑼 𝒊,𝒌 X X 3.79.E+04 2.31.E+03

Once the objective function weights are computed a new NMPC problem is solved minimizing 𝑛

𝑜 a composed objective function of the form ∑𝑖=1 𝑊𝑖 𝜙𝑖 and its optimal solution is used to calculate

the KPIs for the 𝑘 iteration of the Phase II of the proposed algorithm.

References (1)

Garriga, J. L.; Soroush, M. Model Predictive Control Tuning Methods: A Review. Ind. Eng. Chem. Res. 2010, 49 (8), 3505.


41


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(2)

Page 42 of 45

Grune, L.; Jurgen, P. Nonlinear Model Predictive Control. Theory and Algorithms; Isidori, A., Van Schuppen, J. H., Sontag, E. D., Thoma, M., Krstic, M., Eds.; Springer, 2011.

(3)

Tahir, F.; Ohtsuka, T.; Shen, T. Tuning of Nonlinear Model Predictive Controller for the Speed Control of Spark Ignition Engines. Automatic Control Conference (CACS), 2013 CACS International. 2013, pp 216–220.

(4)

Qin, S. J.; Badgwell, T. A. A Survey of Industrial Model Predictive Control Technology. Control Eng. Pract. 2003, 11 (7), 733.

(5)

Shridhar, R.; Cooper, D. J. A Tuning Strategy for Unconstrained Multivariable Model Predictive Control. Ind. Eng. Chem. Res. 1998, 37 (10), 4003.

(6)

Di Cairano, S.; Bemporad, A. Model Predictive Control Tuning by Controller Matching. Automatic Control, IEEE Transactions on. 2010, pp 185–190.

(7)

Francisco, M.; Vega, P. Automatic Tuning of Model Predictive Controllers Based on Multiobjective Optimization. Lat. Am. Appl. Res. 2010, 40, 255.

(8)

Vallerio, M.; Van Impe, J.; Logist, F. Tuning of NMPC Controllers via Multi-Objective Optimisation. Comput. Chem. Eng. 2014, 61 (0), 38.

(9)

Camacho, E. F.; Bordons, C. Model Predictive Control; Grimble, M. J., Johnson, M. A., Eds.; Springer: London, 1999.

(10)

Meadows, E.; Rawlings, J. Model Predictive Control. In Nonlinear Process Control; 1997; pp 233–310.


42

Page 43 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(11)

Susuki, R.; Kawai, F.; Nakazawa, C.; Matsui, T.; Aiyoshi, E. Parameter Optimization of Model Predictive Control Using PSO. SICE Annual Conference, 2008. 2008, pp 1981– 1988.

(12)

Rossiter, J. A. Model-Based Predictive Control. A Practical Approach; Bishop, R. H., Ed.; CRC PRESS: New York, 2005.

(13)

Marlin, T. E. Process Control: Designing Processes and Control Systems for Dynamic Performance, 2nd ed.; McGraw-Hill: Boston, 2000.

(14)

Astrom, K. J.; Wittenmark, B. Computer Controlled Systems. Theory and Design; Pretince Hall Information and System Sciences Series, 1997.

(15)

Botelho, V.; Trierweiler, J. O.; Farenzena, M.; Duraiski, R. Methodology for Detecting Model–Plant Mismatches Affecting Model Predictive Control Performance. Ind. Eng. Chem. Res. 2015, 54 (48), 12072.

(16)

Zavala, V.; Biegler, L. Nonlinear Programming Strategies for State Estimation and Model Predictive Control. In Nonlinear Model Predictive Control; Magni, L., Raimondo, D., Allgöwer, F., Eds.; Lecture Notes in Control and Information Sciences; Springer Berlin Heidelberg; Vol. 384, pp 419–432.

(17)

Jäschke, J.; Yang, X.; Biegler, L. T. Fast Economic Model Predictive Control Based on NLP-Sensitivities. J. Process Control 2014, 24 (8).

(18)

Nafsun, A. I.; Yusoff, N. Effect of Model-Plant Mismatch on MPC Controller Performance. J. Appl. Sci. 2011, 11, 3579.

(19)

Montgomery, D. C. Design and Analysis of Experiments, 7th ed.; John Wiley & Sons, Inc,


43


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 45

2008. (20)

Zavala, V. M. Real-Time Resolution of Conflicting Objectives in Building Energy Management: An Utopia-Tracking Approach. 2012.

(21)

Zavala, V. M.; Flores-Tlacuahuac, A. Stability of Multiobjective Predictive Control: A Utopia-Tracking Approach. Automatica 2012, 48 (10), 2627.

(22)

Flores-Tlacuahuac, A.; Morales, P.; Rivera-Toledo, M. Multiobjective Nonlinear Model Predictive Control of a Class of Chemical Reactors. Ind. Eng. Chem. Res. 2012, 51 (17), 5891.

(23)

Biegler, L. T. Nonlinear Programming. Concepts, Algorithms and Applications to Chemical Process; Optimization, M.-S. S. on, Ed.; SIAM, Society for Industrial and Applied Mathematics: Philadelphi, 2010.

(24)

Bemporad, A.; Muñoz de la Peña, D. Multiobjective Model Predictive Control. Automatica 2009, 45 (12), 2823.

(25)

Pirnay, H.; López-Negrete, R.; Biegler, L. Optimal Sensitivity Based on IPOPT. Math. Program. Comput. 2012, 4 (4), 307.

(26)

Biegler, L. T.; Yang, X.; Fischer, G. A. G. Advances in Sensitivity-Based Nonlinear Model Predictive Control and Dynamic Real-Time Optimization. J. Process Control 2015, 30, 104.

(27)

Huang, R.; Zavala, V. M.; Biegler, L. T. Advanced Step Nonlinear Model Predictive Control for Air Separation Units. J. Process Control 2009, 19 (4), 678.


44

Page 45 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(28)

Biegler, L. T. An Overview of Simultaneous Strategies for Dynamic Optimization. Chem. Eng. Process. Process Intensif. 2007, 46 (11), 1043.

(29)

Lopez-Negrete, R.; D’Amato, F. J.; Biegler, L. T.; Kumar, A. Fast Nonlinear Model Predictive Control: Formulation and Industrial Process Applications. Comput. Chem. Eng. 2013, 51, 55.

(30)

Wächter, A.; Biegler, T. L. On the Implementation of an Interior-Point Filter Line-Search Algorithm for Large-Scale Nonlinear Programming. Math. Program. 2006, 106 (1), 25.

(31)

Hart, W. Python Optimization Modeling Objects (Pyomo). In Operations Research and Cyber-Infrastructure; Chinneck, J., Kristjansson, B., Saltzman, M., Eds.; Operations Research/Computer Science Interfaces; Springer US; Vol. 47, pp 3–19.

(32)

Fogler, H. S. Elements of Chemical Reaction Engineering, 4th ed.; Prentice Hall - PTR: Massachusetts, 2006.

(33)

Santamaría, F. L.; Gómez, J. M. Economic Oriented NMPC for an Extractive Distillation Column Using an Index Hybrid DAE Model Based on Fundamental Principles. Ind. Eng. Chem. Res. 2015, 54 (24), 6344.


45

An Algorithm for Tuning NMPC Controllers with Application to

Recommend Documents