Robust MPC with Output Feedback of Integrating Systems

Abstract

In this work, it is presented a new contribution to the design of a robust MPC with output feedback, input constraints, and uncertain model. Multivariable predictive controllers have been used in industry to reduce the variability of the process output and to allow the operation of the system near to the constraints, where it is usually located the optimum operating point. For this reason, new controllers have been developed with the objective of achieving better performance, simpler control structure, and robustness with respect to model uncertainty. In this work, it is proposed a model predictive controller based on a nonminimal state space model where the state is perfectly known. It is an infinite prediction horizon controller, and it is assumed that there is uncertainty in the stable part of the model, which may also include integrating modes that are frequently present in the process plants. The method is illustrated with a simulation example of the process industry using linear models based on a real process.

1. Introduction

Model predictive control has achieved a remarkable popularity in the process industry with thousands of practical applications [1]. One of the reasons for this industrial acceptance is the ability of MPC to incorporate constraints in the control problem. However, one additional desirable characteristic, still not attended by commercial MPC packages, is closed-loop stability in the presence of model uncertainty, which is usually related to the nonlinearity of the real system. When model uncertainty is considered in the synthesis of the model predictive control, the majority of existing robust algorithms usually demand a computer effort that is prohibitive for practical implementation [2]. From the application viewpoint, we could not find in the control literature a satisfactory solution to the robust MPC problem with output feedback and input constraints. In Rodrigues and Odloak [3], it is presented a formulation to the robust unconstrained MPC with output feedback where stability is achieved through the explicit inclusion of a Lyapunov inequality into the control optimization problem. Later, in Rodrigues and Odloak [4] the method was extended to allow the switching of active input constraints during transient conditions. In Perez [5] the non-minimal state space model (realigned model) proposed in Maciejowski [6] was extended to the incremental form, and this model was used in the controller proposed in Rodrigues and Odloak [4] to produce a robust controller with reduced computational effort by avoiding the use of a state observer, since the state of the realigned model is directly read from the plant. In this approach, the MPC solves the optimization problem that produces the control law in two stages. The first stage is solved offline and computes a family of linear controllers in which the Lyapunov inequalities are added to the control problem so as to ensure state contraction (stability). These linear controllers match all possible configurations of saturation of manipulated variables for a given set of variables that need to be controlled. The second stage is solved online and computes a convex combination of these linear controllers through an optimization problem that includes all the input bounds. These combinations span all possible control configurations in terms of controlled outputs and manipulated inputs. For each control configuration, the linear controller is robust to all process models defining the multiplant uncertainty. The main disadvantage of this approach is that the number of possible control configurations corresponding to the possible combinations of saturated inputs may become very high in systems with large dimension. Then, the off-line step of the control problem may be computationally very expensive or there may not be a feasible solution. The realigned model was applied in González et al. [7] to the development of an infinite horizon MPC nominally stable to open-loop stable and integrating systems. The same non-minimal model was applied in González and Odloak [8] to the development of a robust MPC to uncertain stable systems.

The main objective of this work is to extend the controller proposed in González and Odloak [8] to consider systems with stable and integrating outputs. The method proposed here is based on the non-minimal state space model formulation presented in Maciejowski [6], and the control problem considers that one may have uncertainties in the stable part of the model, while the model for the integrating part is perfectly known. In Section 2, the non-minimal (realigned) model in the incremental form is developed for systems with stable and integrating outputs. The realigned model is split in two parts, one related to the stable outputs and the other related to the integrating outputs. Model uncertainty is only related to the stable outputs, and it is characterized by a family of possible linear time-invariant models. It is assumed that the real model of the system is unknown, but at any time instant the system can be represented by one of the models of this family of models. In practice the models of the integrating outputs are less affected by uncertainty, and here, this uncertainty is neglected. In Section 3, the optimization problem that produces the robust controller is presented. Since the inputs affect both the stable and integrating outputs, the cost function of the robust MPC proposed here includes a combination of the weighted predicted errors in all the outputs. Several results concerning the robustness of the convergence and stability of the resulting closed-loop system are provided. Next, in Section 4, a simulation example, based on linear models of a process system from the oil refining industry, is used to illustrate the application of the new method and finally, in Section 5, the paper is concluded.

2. The Nonminimal State Space Model with Measured State

Consider the following discrete time-invariant model: 𝑦(𝑘)+^𝑛𝑎_𝑖=1𝐴_𝑖𝑦(𝑘−𝑖)=^𝑛𝑏_𝑖=1𝐵_𝑖𝑢(𝑘−𝑖),(1) where 𝑛𝑎 and 𝑛𝑏 are the orders of the discrete model, 𝐴_𝑖∈ℜ^{𝑛𝑦×𝑛𝑦} and 𝐵_𝑖∈ℜ^{𝑛𝑦×𝑛𝑢} are the coefficients that correspond to the parameters of the model, 𝑦(𝑘)∈ℜ^𝑛𝑦 is the system output, and 𝑢(𝑘)∈ℜ^𝑛𝑢 is the system input. Here, it assumed that in the model represented above there are no zero-pole cancellations. In Maciejowski [6], it is shown that this model can be written in a state space form where the state is built with the measurements of the plant inputs and outputs at different time instants. In González et al. [7], it is adopted the same model written in the input incremental form: ⎡⎢⎢⎣𝑥_𝑦(𝑘+1)𝑥_Δ𝑢(𝑘+1)⎤⎥⎥⎦=⎡⎢⎢⎣𝐴_𝑦𝐴_Δ𝑢0_{_}𝐼_{_}⎤⎥⎥⎦⎡⎢⎢⎣𝑥_𝑦(𝑘)𝑥_Δ𝑢(𝑘)⎤⎥⎥⎦+⎡⎢⎢⎣𝐵_Δ𝑢𝐼⎤⎥⎥⎦Δ𝑢(𝑘),𝑦(𝑘)=𝐶_𝑦𝐶_Δ𝑢⎡⎢⎢⎣𝑥_𝑦(𝑘)𝑥_Δ𝑢(𝑘)⎤⎥⎥⎦,(2) where 𝑥_𝑦(𝑘)=𝑦(𝑘)^𝑇𝑦(𝑘−1)^𝑇⋯𝑦(𝑘−𝑛𝑎+1)^𝑇𝑦(𝑘−𝑛𝑎)^𝑇^𝑇∈ℝ^{(𝑛𝑎+1)𝑛𝑦},𝑥_Δ𝑢(𝑘)=Δ𝑢(𝑘−1)^𝑇Δ𝑢(𝑘−2)^𝑇⋯Δ𝑢(𝑘−𝑛𝑏+1)^𝑇^𝑇∈ℝ^{(𝑛𝑏−1)𝑛𝑢},(3)𝑘 is the present sampling instant, Δ𝑢(𝑘)=𝑢(𝑘)−𝑢(𝑘−1) is the input increment, Δ𝑢∈ℜ^𝑛𝑢, 𝑦∈ℜ^𝑛𝑦 is the output. The matrices involved in this model are (see more details in Maciejowski [6])𝐴_𝑦=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝐼−𝐴₁𝐴₁−𝐴₂𝐴₂−𝐴₃𝐴₃−𝐴₄⋯𝐴_𝑛𝑎−1−𝐴_𝑛𝑎𝐴_𝑛𝑎𝐼_{𝑛𝑦×𝑛𝑦}000⋯000𝐼_{𝑛𝑦×𝑛𝑦}00⋯0000𝐼_{𝑛𝑦×𝑛𝑦}0⋯00000𝐼_{𝑛𝑦×𝑛𝑦}⋯00⋮⋮⋮⋮⋱⋮⋮0000⋯𝐼_{𝑛𝑦×𝑛𝑦}0⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,𝐴_𝑦∈ℝ^{(𝑛𝑎+1)𝑛𝑦×(𝑛𝑎+1)𝑛𝑦},𝐴_Δ𝑢=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣𝐵₂⋯𝐵_𝑛𝑏0_{𝑛𝑦×𝑛𝑢}⋮0_{𝑛𝑦×𝑛𝑢}0_{𝑛𝑦×𝑛𝑢}⋮0_{𝑛𝑦×𝑛𝑢}⋮⋱⋮0_{𝑛𝑦×𝑛𝑢}⋯0_{𝑛𝑦×𝑛𝑢}⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦,𝐴_Δ𝑢∈ℝ^{(𝑛𝑎+1)⋅𝑛𝑦×(𝑛𝑏−1)⋅𝑛𝑢},𝐼_{_}=⎡⎢⎢⎢⎢⎢⎢⎣0_𝑛𝑢⋯0_𝑛𝑢𝐼_𝑛𝑢⋯0_𝑛𝑢⋮⋱⋮0_𝑛𝑢𝐼_𝑛𝑢0_𝑛𝑢⎤⎥⎥⎥⎥⎥⎥⎦,𝐼_{_}∈ℝ^{(𝑛𝑏−1)𝑛𝑢×(𝑛𝑏−1)𝑛𝑢},0_{_}∈ℝ^{(𝑛𝑏−1)𝑛𝑢×(𝑛𝑎+1)𝑛𝑦},𝐵_Δ𝑢=⎡⎢⎢⎣𝐵₁0⎤⎥⎥⎦=⎡⎢⎢⎢⎢⎢⎢⎣𝐵₁0_{𝑛𝑦×𝑛𝑢}⋮0_{𝑛𝑦×𝑛𝑢}⎤⎥⎥⎥⎥⎥⎥⎦,𝐵₁∈ℝ^{𝑛𝑦×𝑛𝑢},𝐵_Δ𝑢∈ℝ^{[(𝑛𝑎+1)⋅𝑛𝑦]×𝑛𝑢},𝐼=⎡⎢⎢⎢⎢⎢⎢⎣𝐼_𝑛𝑢0_𝑛𝑢⋮0_𝑛𝑢⎤⎥⎥⎥⎥⎥⎥⎦,𝐼∈ℝ^{[(𝑛𝑏−1)⋅𝑛𝑢]×𝑛𝑢},𝐶_𝑦=⎡⎢⎢⎢⎣𝐼_{𝑛𝑦×𝑛𝑦}0_{𝑛𝑦×𝑛𝑦}⋯0_{𝑛𝑦×𝑛𝑦𝑛𝑎+1}⎤⎥⎥⎥⎦,𝐶_Δ𝑢=⎡⎢⎢⎢⎣0_{𝑛𝑦×𝑛𝑢}0_{𝑛𝑦×𝑛𝑢}⋯0_{𝑛𝑦×𝑛𝑢𝑛𝑏−1}⎤⎥⎥⎥⎦.(4)In the system defined by (2), the input and output values read from the plant are used to realign the model, which gives a better disturbance rejection capability to the controller. In González and Odloak [8], it is shown that the model presented above is detectable and stabilizable. More detailed discussions of minimal state space representations and the resulting observability and controllability of time-delayed systems can be found in De La Sen [9].

For control implementation, in order to better locate the uncertainties along the process model, it is adequate to divide the model in two separated parts. The first part of the model is related to the pure integrating outputs, and the second part is related to the stable outputs. Here, it is assumed that the process has no outputs related simultaneously with stable and integrating modes. Then, suppose that the model defined by (2) is written for the integrating outputs: ⎡⎢⎢⎣𝑥^𝑖_𝑦(𝑘+1)𝑥^𝑖_Δ𝑢(𝑘+1)⎤⎥⎥_{𝑥𝑖(𝑘+1)}⎦=⎡⎢⎢⎣𝐴^𝑖_𝑦𝐴^𝑖_Δ𝑢0_{_}𝐼_{_}⎤⎥⎥_{𝐴𝑖}⎦⎡⎢⎢⎣𝑥^𝑖_𝑦(𝑘)𝑥^𝑖_Δ𝑢(𝑘)⎤⎥⎥_{𝑥𝑖(𝑘)}⎦+⎡⎢⎢⎣𝐵^𝑖_Δ𝑢𝐼⎤⎥⎥_{𝐵𝑖}⎦Δ𝑢(𝑘),(5)𝑦^𝑖(𝑘)=𝐶^𝑖_𝑦𝐶^𝑖_Δ𝑢⎡⎢⎢⎣𝑥^𝑖_𝑦(𝑘)𝑥^𝑖_Δ𝑢(𝑘)⎤⎥⎥⎦.(6) Analogously, a similar model can be written for the stable outputs: ⎡⎢⎢⎣𝑥^𝑠_{𝑦,𝜃𝑝}(𝑘+1)𝑥^𝑠_Δ𝑢(𝑘+1)⎤⎥⎥_{𝑥𝑠𝜃𝑝(𝑘+1)}⎦=⎡⎢⎢⎣𝐴^𝑠_𝑦𝜃_𝑝𝐴^𝑠_Δ𝑢𝜃_𝑝0_{_}𝐼_{_}⎤⎥⎥_{𝐴𝑠𝜃𝑝}⎦⎡⎢⎢⎣𝑥^𝑠_{𝑦,𝜃𝑝}(𝑘)𝑥^𝑠_Δ𝑢(𝑘)⎤⎥⎥_{𝑥𝑠𝜃𝑝(𝑘)}⎦+⎡⎢⎢⎣𝐵^𝑠_Δ𝑢𝜃_𝑝𝐼⎤⎥⎥_{𝐵𝑠𝜃𝑝}⎦Δ𝑢(𝑘),(7)𝑦^𝑠_𝜃𝑝(𝑘)=𝐶^𝑠_𝑦𝐶^𝑠_Δ𝑢⎡⎢⎢⎣𝑥^𝑠_{𝑦,𝜃𝑝}(𝑘)𝑥^𝑠_Δ𝑢(𝑘)⎤⎥⎥⎦.(8) Observe that some of the components of the model defined in (7) and (8) depend on 𝜃_𝑝, that is, the vector of parameters that define the model (𝜃_𝑝=[𝐴_𝑝,𝐵_𝑝]). These parameters are unknown, which characterizes the uncertainty in the model related to the stable outputs. The simplest form to consider uncertainty is through the multimodel uncertainty where 𝜃_𝑝 belongs to a finite set 𝜃_𝑝∈Ω={𝜃₁,…,𝜃_𝐿}. This means that the true model of the process is unknown, but it is one of the 𝐿 models of a finite set of known models.

Within the conventional MPC formulation, the control horizon 𝑚 is such that Δ𝑢(𝑘+𝑗)=0, 𝑗=𝑚,𝑚+1,…, then, with the model representation defined in (2) and (3), it is easy to show that 𝑥_Δ𝑢(𝑘+𝑛+𝑗)=0,𝑗=1,2,…, where 𝑛=𝑚+𝑛𝑏−1 and 𝑥_𝑦(𝑘+𝑛+𝑗)=𝐴_𝑦𝑥_𝑦(𝑘+𝑛+𝑗−1),𝑗=1,2,….

Then, after time step 𝑘+𝑛, it is possible to apply the Jordan matrix decomposition of 𝐴_𝑦 as follows: 𝐴_𝑦𝑉=𝑉𝐴_𝑑,(9) where 𝐴_𝑑 is a block diagonal matrix where the eigenvalues of 𝐴_𝑦 appear in its main diagonal. If 𝐴_𝑦 is a full-rank square matrix, then 𝑉 is invertible and it is possible to make a change of variables such that the integrating states related to the incremental form of the model are separated from the remaining states of the system. For instance, for the integrating system defined in (5) and (6), one can define 𝑥^𝑖_𝑦(𝑘)=𝑉_𝑖𝑧^𝑖(𝑘),𝑥^𝑖_𝑦(𝑘)=𝑉^Δ_𝑖𝑉^𝑖_𝑖⎡⎢⎢⎣𝑧^Δ_𝑖(𝑘)𝑧^𝑖_𝑖(𝑘)⎤⎥⎥⎦,(10) where 𝑧^Δ_𝑖 corresponds to the integrating states related to the incremental form of the model and 𝑧^𝑖_𝑖 is the state associated with the real integrating states of the system. 𝑉_𝑖 is obtained from the Jordan decomposition of matrix 𝐴^𝑖_𝑦 as follows: 𝐴^𝑖_𝑦𝑉_𝑖=𝑉_𝑖𝐴^𝑖_𝑑.(11) It can be shown that 𝐴^𝑖_𝑑=^{𝐼𝑛𝑦𝑖𝐼𝑛𝑦𝑖0𝐼}_𝑛𝑦𝑖 where 𝑛𝑦𝑖 is the number of integrating outputs.

It should be noted that, after time step 𝑘+𝑛, the transformed state of the integrating system will evolve according to the equation ⎡⎢⎢⎣𝑧^Δ_𝑖(𝑘+𝑛+1)𝑧^𝑖_𝑖(𝑘+𝑛+1)⎤⎥⎥⎦=⎡⎢⎢⎣𝐼_𝑛𝑦𝐼_𝑛𝑦0𝐼_𝑛𝑦⎤⎥⎥⎦⎡⎢⎢⎣𝑧^Δ_𝑖(𝑘+𝑛)𝑧^𝑖_𝑖(𝑘+𝑛)⎤⎥⎥⎦.(12) Analogously, for the system represented through (7) and (8), the following state transformation can be defined: 𝑥^𝑠_{𝑦,𝜃𝑝}(𝑘)=𝑉_{𝑠,𝜃𝑝}𝑧^𝑠_𝜃𝑝(𝑘)=𝑉^Δ_{𝑠,𝜃𝑝}𝑉^𝑠_{𝑠,𝜃𝑝}⎡⎢⎢⎣𝑧^Δ_{𝑠,𝜃𝑝}(𝑘)𝑧^𝑠_{𝑠,𝜃𝑝}(𝑘)⎤⎥⎥⎦,(13) where 𝑧^Δ_{𝑠,𝜃𝑝} corresponds to the integrating states related to the incremental form of the model and 𝑧^𝑠_{𝑠,𝜃𝑝} is the state associated with the stable states of the system corresponding to the parameters represented as 𝜃_𝑝. 𝑉_{𝑠,𝜃𝑝} is obtained from the Jordan decomposition of matrix 𝐴^𝑠_𝑦(𝜃_𝑝) as follows: 𝐴^𝑠_𝑦𝜃_𝑝𝑉_{𝑠,𝜃𝑝}=𝑉_{𝑠,𝜃𝑝}𝐴^𝑠_{𝑑,𝜃𝑝}.(14) It can be shown that 𝐴^𝑠_{𝑑,𝜃𝑝}=⎡⎢⎢⎣𝐼_𝑛𝑦𝑠00𝐹^𝑠_𝜃𝑝⎤⎥⎥⎦(15) where 𝑛𝑦𝑠 is the number of stable outputs of the system. If the stable poles of the system are nonrepeated, 𝐹^𝑠_𝜃𝑝 is a diagonal matrix containing the eigenvalues of 𝐴^𝑠_𝑦(𝜃_𝑝).

After time step 𝑘+𝑛, the transformed state corresponding to the stable outputs will evolve according to the equation ⎡⎢⎢⎣𝑧^Δ_{𝑠,𝜃𝑝}(𝑘+𝑛+1)𝑧^𝑠_{𝑠,𝜃𝑝}(𝑘+𝑛+1)⎤⎥⎥⎦=⎡⎢⎢⎣𝐼_𝑛𝑦𝑠00𝐹^𝑠_𝜃𝑝⎤⎥⎥⎦⎡⎢⎢⎣𝑧^Δ_{𝑠,𝜃𝑝}(𝑘+𝑛)𝑧^𝑠_{𝑠,𝜃𝑝}(𝑘+𝑛)⎤⎥⎥⎦.(16)

3. The Robust MPC with Output Feedback

The system outputs need to be controlled through the manipulation of the nu inputs. Then, the robust MPC proposed here is based on the following objective function: 𝐽_𝑘𝜃_𝑝=^𝑛_𝑗=1𝑦^𝑖(𝑘+𝑗∣𝑘)−𝑦^𝑠𝑝_𝑖−𝛿^𝑖_Δ,𝑘−𝑗𝛿^𝑖_𝑖,𝑘^𝑇×𝑄_𝑖𝑦^𝑖(𝑘+𝑗∣𝑘)−𝑦^𝑠𝑝_𝑖−𝛿^𝑖_Δ,𝑘−𝑗𝛿^𝑖_𝑖,𝑘+^∞_𝑗=𝑛+1𝑦^𝑖(𝑘+𝑗∣𝑘)−𝑦^𝑠𝑝_𝑖−𝛿^𝑖_Δ,𝑘−𝑗𝛿^𝑖_𝑖,𝑘^𝑇×𝑄_𝑖𝑦^𝑖(𝑘+𝑗∣𝑘)−𝑦^𝑠𝑝_𝑖−𝛿^𝑖_Δ,𝑘−𝑗𝛿^𝑖_𝑖,𝑘+^𝑛_𝑗=1𝑦^𝑠_𝜃𝑝(𝑘+𝑗∣𝑘)−𝑦^𝑠𝑝_𝑠−𝛿^𝑠_Δ,𝑘𝜃_𝑝^𝑇×𝑄_𝑠𝑦^𝑠_𝜃𝑝(𝑘+𝑗∣𝑘)−𝑦^𝑠𝑝_𝑠−𝛿^𝑠_Δ,𝑘𝜃_𝑝+^∞_𝑗=𝑛+1𝑦^𝑠_𝜃𝑝(𝑘+𝑗∣𝑘)−𝑦^𝑠𝑝_𝑠−𝛿^𝑠_Δ,𝑘𝜃_𝑝^𝑇×𝑄_𝑠𝑦^𝑠_𝜃𝑝(𝑘+𝑗∣𝑘)−𝑦^𝑠𝑝_𝑠−𝛿^𝑠_Δ,𝑘𝜃_𝑝+^𝑚−1_𝑗=0Δ𝑢(𝑘+𝑗∣𝑘)^𝑇𝑅Δ𝑢(𝑘+𝑗∣𝑘)+𝛿^𝑖_Δ,𝑘^𝑇𝑆^𝑖_Δ𝛿^𝑖_Δ,𝑘+𝛿^𝑠_Δ,𝑘𝜃_𝑝^𝑇𝑆^𝑠_Δ𝛿^𝑠_Δ,𝑘𝜃_𝑝+𝛿^𝑖_𝑖,𝑘^𝑇𝑆^𝑖_𝑖𝛿^𝑖_𝑖,𝑘,(17) where 𝑦^𝑖(𝑘+𝑗∣𝑘) is prediction of integrating output at time 𝑘+𝑗 computed at time 𝑘 and 𝑦^𝑠_𝜃𝑝(𝑘+𝑗∣𝑘) is the prediction of the stable output corresponding to model 𝜃_𝑝; 𝛿^𝑖_Δ,𝑘, 𝛿^𝑠_Δ,𝑘, and 𝛿^𝑖_𝑖,𝑘 are slack variables that are included in the control problem to guarantee that 𝐽_𝑘(𝜃_𝑝) is bounded. 𝑄_𝑖, 𝑄_𝑠, 𝑅, 𝑆^𝑖_Δ, 𝑆^𝑠_Δ and 𝑆^𝑖_𝑖 are positive weighting matrices. In the second term on the right-hand side of (17), the prediction of the integrating output can be written as follows: 𝑦^𝑖(𝑘+𝑛+𝑗∣𝑘)=𝐶^𝑖_𝑦𝑥^𝑖_𝑦(𝑘+𝑛+𝑗∣𝑘)+𝐶^𝑖_Δ𝑢𝑥^𝑖_Δ𝑢(𝑘+𝑛+𝑗∣𝑘).(18) But, since 𝑥^𝑖_Δ𝑢(𝑘+𝑛+𝑗∣𝑘)=0, (18) becomes 𝑦^𝑖(𝑘+𝑛+𝑗∣𝑘)=𝐶^𝑖_𝑦𝑥^𝑖_𝑦(𝑘+𝑛+𝑗∣𝑘).(19) Then, using (10) and (12), (19) can be written as follows: 𝑦^𝑖(𝑘+𝑛+𝑗∣𝑘)=𝐶^𝑖_𝑦𝑉^Δ_𝑖𝑧^Δ_𝑖(𝑘+𝑛∣𝑘)+𝑗𝑧^𝑖_𝑖(𝑘+𝑛∣𝑘)+𝐶^𝑖_𝑦𝑉^𝑖_𝑖𝑧^𝑖_𝑖(𝑘+𝑛∣𝑘).(20) If (20) is substituted into the infinite sum represented by the second term on the right-hand side of (17), it is easy to see that the objective function defined in (17) will be bounded only if the following constraints are satisfied: 𝐶^𝑖_𝑦𝑉^Δ_𝑖𝑧^Δ_𝑖(𝑘+𝑛∣𝑘)+𝐶^𝑖_𝑦𝑉^𝑖_𝑖𝑧^𝑖_𝑖(𝑘+𝑛∣𝑘)−𝑦^𝑠𝑝_𝑖−𝛿^𝑖_Δ,𝑘−𝑛𝛿^𝑖_𝑖,𝑘=0,(21)𝐶^𝑖_𝑦𝑉^Δ_𝑖𝑧^𝑖_𝑖(𝑘+𝑛∣𝑘)−𝛿^𝑖_𝑖,𝑘=0.(22) Analogously, in the infinite sum represented by the fourth term on the right-hand side of (17), the prediction of the stable output can be written as follows: 𝑦^𝑠_𝜃𝑝(𝑘+𝑛+𝑗∣𝑘)=𝐶^𝑠_𝑦𝑥^𝑠_{𝑦,𝜃𝑝}(𝑘+𝑛+𝑗∣𝑘).(23) Now, using (13) and (16), equation (23) can be written as follows: 𝑦^𝑠_𝜃𝑝(𝑘+𝑛+𝑗∣𝑘)=𝐶^𝑠_𝑦𝑉^Δ_{𝑠,𝜃𝑝}𝑧^Δ_{𝑦,𝜃𝑝}(𝑘+𝑛∣𝑘)+𝐶^𝑠_𝑦𝑉^𝑠_{𝑠,𝜃𝑝}𝐹^𝑠_𝜃𝑝^𝑗𝑧^𝑠_{𝑦,𝜃𝑝}(𝑘+𝑛∣𝑘).(24) If the expression for the output prediction obtained in (24) is substituted into the infinite sum represented by the fourth term on the right-hand side of (17), it is easy to see that the cost function will be bounded only if the following constraint is satisfied: 𝐶^𝑠_𝑦𝑉^Δ_{𝑠,𝜃𝑝}𝑧^Δ_{𝑦,𝜃𝑝}(𝑘+𝑛∣𝑘)−𝑦^𝑠𝑝_𝑠−𝛿^𝑠_Δ,𝑘𝜃_𝑝=0.(25) Consequently, the optimization problem that produces the control law proposed here will have to include the constraints defined in (21), (22), and (25).

Equation (21) can be developed in terms of the available states at time 𝑘 and the vector of the future control moves as follows: 𝐶^𝑖_𝑦𝑉^Δ_𝑖𝑁^𝑖₁+𝐶^𝑖_𝑦𝑉^𝑖_𝑖𝑁^𝑖₂𝑉⁻¹_𝑖𝑁^𝑖𝑥^𝑖(𝑘+𝑛∣𝑘)−𝑦^𝑠𝑝_𝑖−𝛿^𝑖_Δ,𝑘−𝑛𝛿^𝑖_𝑖,𝑘=0,(26) where 𝑁^𝑖₁ is a matrix of ones and zeros that collects component 𝑧^Δ_𝑖 from 𝑧_𝑖 and matrix 𝑁^𝑖₂ collects component 𝑧^𝑖_𝑖 from 𝑧_𝑖. The model defined in (5) can be used to represent (26) in terms of the integrating state that at time 𝑘, is measured as follows: 𝐶^𝑖_𝑦𝑉^Δ_𝑖𝑁^𝑖₁+𝐶^𝑖_𝑦𝑉^𝑖_𝑖𝑁^𝑖₂𝑉⁻¹_𝑖𝑁^𝑖𝐴^𝑖^𝑛−𝑚𝑥^𝑖(𝑘)+𝐵^𝑖_𝑚Δ𝑢_𝑘−𝑦^𝑠𝑝_𝑖−𝛿^𝑖_Δ,𝑘−𝑛𝛿^𝑖_𝑖,𝑘=0,(27) where Δ𝑢_𝑘=[Δ𝑢(𝑘∣𝑘)^𝑇⋯Δ𝑢(𝑘+𝑚−1∣𝑘)^𝑇]^𝑇 and 𝐵^𝑖_𝑚=[(𝐴^𝑖)^𝑛−1𝐵^𝑖⋯(𝐴^𝑖)^𝑛−𝑚𝐵^𝑖].

The constraint (22) can also be represented in terms of the present state and the vector of future control moves: 𝐶^𝑖_𝑦𝑉^Δ_𝑖𝑁^𝑖₁𝑉⁻¹_𝑖𝑁^𝑖𝐴^𝑖^𝑛−𝑚𝑥^𝑖(𝑘)+𝐵^𝑖_𝑚Δ𝑢_𝑘−𝛿^𝑖_𝑖,𝑘=0.(28) Analogously, (25) can be written as follows: 𝐶^𝑠_𝑦𝑉^Δ_{𝑠,𝜃𝑝}𝑁^𝑠₁𝑉⁻¹_{𝑠,𝜃𝑝}𝑁^𝑠₂𝐴^𝑠𝜃_𝑝^𝑛−𝑚𝑥^𝑠_𝜃𝑝(𝑘)+𝐵^𝑠_𝑚𝜃_𝑝Δ𝑢_𝑘−𝑦^𝑠𝑝_𝑠−𝛿^𝑠_Δ,𝑘𝜃_𝑝=0,(29) where 𝑁^𝑠₁ captures the state component 𝑧^Δ_{𝑠,𝜃𝑝} of state 𝑧^𝑠_𝜃𝑝 and 𝑁^𝑠₂ captures the component 𝑥^𝑠_{𝑦,𝜃𝑝} of state 𝑥^𝑠_𝜃𝑝. With the constraints discussed above, the control cost function defined in (17) can be written as follows: 𝐽_𝑘𝜃_𝑝=^𝑛_𝑗=1𝑦^𝑖(𝑘+𝑗∣𝑘)−𝑦^𝑠𝑝_𝑖−𝛿^𝑖_Δ,𝑘−𝑗𝛿^𝑖_𝑖,𝑘^𝑇×𝑄_𝑖𝑦^𝑖(𝑘+𝑗∣𝑘)−𝑦^𝑠𝑝_𝑖−𝛿^𝑖_Δ,𝑘−𝑗𝛿^𝑖_𝑖,𝑘+^𝑛_𝑗=1𝑦^𝑠_𝜃𝑝(𝑘+𝑗∣𝑘)−𝑦^𝑠𝑝_𝑠−𝛿^𝑠_Δ,𝑘𝜃_𝑝^𝑇×𝑄_𝑠𝑦^𝑠_𝜃𝑝(𝑘+𝑗∣𝑘)−𝑦^𝑠𝑝_𝑠−𝛿^𝑠_Δ,𝑘𝜃_𝑝+𝑥^𝑠_𝜃𝑝(𝑘+𝑛∣𝑘)^𝑇(𝑁^𝑠)^𝑇𝑉⁻¹_{𝑠,𝜃𝑝}^𝑇𝑁^𝑠₃^𝑇×𝑄⁻_𝑠𝑁^𝑠₃𝑉⁻¹_{𝑠,𝜃𝑝}𝑁^𝑠𝑥^𝑠_𝜃𝑝(𝑘+𝑛∣𝑘)+^𝑚−1_𝑗=0Δ𝑢(𝑘+𝑗∣𝑘)^𝑇𝑅Δ𝑢(𝑘+𝑗∣𝑘)+𝛿^𝑖_Δ,𝑘^𝑇𝑆^𝑖_Δ𝛿^𝑖_Δ,𝑘+𝛿^𝑠_Δ,𝑘𝜃_𝑝^𝑇𝑆^𝑠_Δ𝛿^𝑠_Δ,𝑘𝜃_𝑝+𝛿^𝑖_𝑖,𝑘^𝑇𝑆^𝑖_𝑖𝛿^𝑖_𝑖,𝑘,(30) where 𝑄⁻_𝑠 is obtained from the solution to the following equation: 𝑄⁻_𝑠−𝐹^𝑠_𝜃𝑝𝑄^𝑇−_𝑠𝐹^𝑠_𝜃𝑝=𝐹^𝑠_𝜃𝑝^𝑇𝑉^𝑠_{𝑠,𝜃𝑝}^𝑇𝐶^𝑠_𝑦^𝑇𝑄_𝑠𝐶^𝑠_𝑦𝑉^𝑠_{𝑠,𝜃𝑝}𝐹^𝑠_𝜃𝑝.(31) Observe that (30) is a finite horizon cost function with a terminal weight computed through (31).

One can now define the robust MPC with output feedback for systems with integrating and stable outputs. The controller is robust in a sense that it maintains stability even in the presence of uncertainty in the model related to the stable outputs. Assuming that 𝜃_𝑁 corresponds to the nominal or most probable model, the proposed controller is based on the solution to the following problem.

Problem 1. Consider the following objective function: min_{Δ𝑢𝑘,𝛿𝑖𝑖,𝑘,𝛿𝑖Δ,𝑘,𝛿𝑠Δ,𝑘(𝜃𝑝)𝑝=1,...,𝐿}𝐽_𝑘𝜃_𝑁(32) subject to (26), (28) and 𝐶^𝑠_𝑦𝑉^Δ_{𝑠,𝜃𝑝}𝑁^𝑠₁𝑉⁻¹_{𝑠,𝜃𝑝}𝑁^𝑠₂𝐴^𝑠𝜃_𝑝^𝑛−𝑚𝑥^𝑠_𝜃𝑝(𝑘)+𝐵^𝑠_𝑚𝜃_𝑝Δ𝑢_𝑘−𝑦^𝑠𝑝_𝑠−𝛿^𝑠_Δ,𝑘𝜃_𝑝=0,𝑝=1,…,𝐿−Δ𝑢_max≤Δ𝑢(𝑘+𝑗∣𝑘)≤Δ𝑢_max,𝑗=0,…,𝑚−1𝑢_min≤𝑢(𝑘+𝑗∣𝑘)≤𝑢_max,𝑗=0,…,𝑚−1𝐽_𝑘𝜃_𝑝≤𝐽_𝑘̃‌𝜃_𝑝,𝑝=1,…,𝐿,(33)‖‖𝛿^𝑖_𝑖,𝑘‖‖≤‖‖𝛿̃‌^𝑖_𝑖,𝑘‖‖,(34) where 𝐽_𝑘̃‌(𝜃_𝑝) is computed with Δ𝑢_𝑘̃=Δ𝑢^∗(𝑘∣𝑘−1)Δ𝑢^∗(𝑘+1∣𝑘−1)⋯Δ𝑢^∗(𝑘+𝑚−2∣𝑘−1)0,(35)where Δ𝑢^∗(𝑘+𝑗∣𝑘−1) corresponds to the optimal solution to Problem 1 at time step 𝑘−1 and 𝛿̃‌^𝑖_𝑖,𝑘,𝛿̃‌^𝑖_Δ,𝑘,𝛿̃‌^𝑠_Δ,𝑘(𝜃_𝑝),𝑝=1,…,𝐿, which are obtained from the solution to the following equations: 𝐶^𝑠_𝑦𝑉^Δ_{𝑠,𝜃𝑝}𝑁^𝑠₁𝑉⁻¹_{𝑠,𝜃𝑝}𝑁^𝑠₂𝐴^𝑠𝜃_𝑝^𝑛−𝑚𝑥^𝑠_𝜃𝑝(𝑘)+𝐵^𝑠_𝑚𝜃_𝑝Δ𝑢_𝑘̃−𝑦^𝑠𝑝_𝑠−𝛿̃‌^𝑠_Δ,𝑘𝜃_𝑝=0,𝑝=1,…,𝐿,𝐶^𝑖_𝑦𝑉^Δ_𝑖𝑁^𝑖₁+𝐶^𝑖_𝑦𝑉^𝑖_𝑖𝑁^𝑖₂𝑉⁻¹_𝑖𝑁^𝑖𝐴^𝑖^𝑛−𝑚𝑥^𝑖(𝑘)+𝐵^𝑖_𝑚Δ𝑢_𝑘̃−𝑦^𝑠𝑝_𝑖−𝛿̃‌^𝑖_Δ,𝑘−𝑛𝛿̃‌^𝑖_𝑖,𝑘=0,𝐶^𝑖_𝑦𝑉^Δ_𝑖𝑁^𝑖₁𝑉⁻¹_𝑖𝑁^𝑖𝐴^𝑖^𝑛−𝑚𝑥^𝑖(𝑘)+𝐵^𝑖_𝑚Δ𝑢_𝑘̃−𝛿̃‌^𝑖_𝑖,𝑘=0.(36)

Remark 2. The constraint defined in (33) corresponds to the conventional cost contracting constraint first proposed in Badgwell [10] in the development of a robust MPC for the regulator problem and extended in Odloak [11] for the output tracking problem. In these works, it is shown that, for open-loop stable systems with multiplant uncertainty, the inclusion of constraint (33) forces the cost function of the true plant to be a Lyapunov function for the closed-loop system.

Remark 3. The constraint defined in (34) is intended to force slack 𝛿^𝑖_𝑖,𝑘 to converge to zero. It will be shown latter that if Problem 1 is feasible with this slack kept equal to zero, the control cost function of the true plant becomes a Lyapunov function for the uncertain system and the controller resulting from the solution to Problem 1. However, if there is uncertainty on the model related to the integrating part of the system, which means that there is uncertainty in matrices 𝐴^𝑖 and/or 𝐵^𝑖, then, it is easy to show that it is not possible to find a solution to Problem 1 that satisfies constraint (28) for the different matrices 𝐴^𝑖 and 𝐵^𝑖 with 𝛿^𝑖_𝑖,𝑘(𝜃_𝑝)=0,𝑝=1,…,𝐿. This means that the approach followed here cannot be applied for systems with uncertainty in the gains of the integrating outputs. In the next theorems we prove the robust stability of the controller obtained through the solution to Problem 1.

The proof of the stability of an MPC usually involves two ingredients: recursive feasibility of the control problem and the existence of a Lyapnvov function for the closed-loop system. The theorem below shows the recursive feasibility of Problem 1, and the other two theorems show that the cost function can be interpreted as a Lyapunov function for the system.

Theorem 4. If Problem 1 has a feasible solution at time step 𝑘, it will remain feasible at any subsequent time step 𝑘+𝑗.

Proof. Suppose that at time 𝑘 the optimum solution is represented as follows: Δ𝑢^∗_𝑘,𝛿^∗𝑖_𝑖,𝑘,𝛿^∗𝑖_Δ,𝑘,𝛿^∗𝑠_Δ,𝑘𝜃_𝑝,𝑝=1,…,𝐿.(37) Then, it is easy to show that, at time step 𝑘+1, the following solution, inherited from the optimum solution at time 𝑘, will be feasible:Δ𝑢_𝑘+1=Δ𝑢^∗(𝑘+1∣𝑘)Δ𝑢^∗(𝑘+2∣𝑘1)⋯Δ𝑢^∗(𝑘+𝑚−1∣𝑘)0,𝛿^𝑖_𝑖,𝑘+1=𝛿̃‌^𝑖_𝑖,𝑘+1,𝛿^𝑖_Δ,𝑘+1=𝛿̃‌^𝑖_Δ,𝑘,𝛿^𝑠_Δ,𝑘+1𝜃_𝑝=𝛿̃‌^𝑠_Δ,𝑘𝜃_𝑝,𝑝=1,…,𝐿,(38)where 𝛿̃‌^𝑖_𝑖,𝑘+1,𝛿̃‌^𝑖_Δ,𝑘+1,𝛿̃‌^𝑠_Δ,𝑘+1(𝜃_𝑝),𝑝=1,...,𝐿 are computed as defined in the enunciate of the theorem and consequently, one can always find a solution to Problem 1.

Since Problem 1 has recursive feasibility, the next theorem shows that the consideration of constraint (34) in the control problem guarantees that, if the system is not disturbed, 𝛿^𝑖_𝑖,𝑘 that is the slack related to the integrating outputs will converge to zero.

Theorem 5. If the system defined in (5), (7), (6), and (8) is stabilizable, the solution to Problem 1 at successive time steps drives the slacks 𝛿^𝑖_𝑖,𝑘 related to the integrating outputs to zero.

Proof. Suppose that, at time step 𝑘, Problem 1 is solved and the solution defined in (37) is obtained. Then, since it is assumed that there is no uncertainty in the model related to the integrating outputs, it is easy to see that for the undisturbed system 𝛿̃‌^𝑖_Δ,𝑘+1=𝛿^∗𝑖_Δ,𝑘, and consequently, ‖𝛿^∗𝑖_Δ,𝑘+1‖≤‖𝛿^∗𝑖_Δ,𝑘‖. Then, if the system is stabilizable, 𝛿^𝑖_Δ,𝑘 will be driven to zero.

Theorem 6. Consider the system represented in (5), (6), (7), and (8). If the conditions defined in Theorems 4 and 5 are satisfied and the desired steady state is reachable, then, the solution of Problem 1 at successive time steps will produce a sequence of control moves that will drive the system outputs to the desired targets.

Proof. Suppose that the convergence of 𝛿^𝑖_𝑖,𝑘 to zero has already been achieved and Problem 1 is solved at time step 𝑘. Let us designate the resulting optimum solution as Δ𝑢^∗_𝑘,𝛿^∗𝑖_𝑖,𝑘=0,𝛿^∗𝑖_Δ,𝑘,𝛿^∗𝑠_Δ,𝑘𝜃_𝑝,𝑝=1,…,𝐿.(39) Let us now define the set of parameters corresponding to the true model as 𝜃_𝑇. The resulting optimum control move Δ𝑢^∗(𝑘∣𝑘) is implemented in the true system and we move to time instant 𝑘+1 where Problem 1 is solved again. At this time step, as shown in Theorem 4, Δ𝑢̃_𝑘+1,𝛿̃‌^𝑖_𝑖,𝑘+1,𝛿̃‌^𝑖_Δ,𝑘+1,𝛿̃‌^𝑠_Δ,𝑘+1(𝜃_𝑝),𝑝=1,…,𝐿 is a feasible solution to Problem 1. Also, since 𝛿̃‌^𝑖_𝑖,𝑘+1=𝛿^∗𝑖_𝑖,𝑘=0, it is easy to show that 𝛿̃‌^𝑖_Δ,𝑘+1=𝛿^∗𝑖_Δ,𝑘,𝛿̃‌^𝑠_Δ,𝑘+1(𝜃_𝑇)=𝛿^∗𝑠_Δ,𝑘(𝜃_𝑇). Thus, it is easy to see that 𝐽̃‌_𝑘+1(𝜃_𝑇)=𝐽^∗_𝑘(𝜃_𝑇), and consequently, 𝐽^∗_𝑘+1(𝜃_𝑇)≤𝐽^∗_𝑘(𝜃_𝑇). Now, since 𝐽_𝑘(𝜃_𝑇) is positive and bounded below by zero, it can be interpreted as a Lyapunov function for the closed-loop system with the proposed controller. This means that 𝑦^𝑖(𝑘+𝑗∣𝑘)−𝑦^𝑠𝑝_𝑖−𝛿^𝑖_Δ,𝑘 and 𝑦^𝑠_𝜃𝑇(𝑘+𝑗∣𝑘)−𝑦^𝑠𝑝_𝑠−𝛿^𝑠_Δ,𝑘(𝜃_𝑇) will converge to zero as 𝑗→∞. Following the same steps as in González et al. [12], it can be shown that if weights 𝑆^𝑖_Δ and 𝑆^𝑠_Δ are large enough in comparison to 𝑅 and the desired set point is reachable, slacks 𝛿^𝑖_Δ,𝑘 and 𝛿^𝑠_Δ,𝑘(𝜃_𝑇) will converge to zero and the outputs will converge to the set points. If the targets are not reachable, the distance to the targets will be minimized according to the norm (𝛿^𝑖_Δ,𝑘)^𝑇𝑆^𝑖_Δ𝛿^𝑖_Δ,𝑘+(𝛿^𝑠_Δ,𝑘(𝜃_𝑇))^𝑇𝑆^𝑠_Δ𝛿^𝑠_Δ,𝑘(𝜃_𝑇).

4. Simulation Example: The Distillation Column System

The proposed control strategy was tested with a small dimension system of the process industry. The system is part of a distillation column where isobutene is separated from n-butane in the alkylation unit of an oil refinery. This system was studied in Carrapiço et al. [13] that implemented an infinite horizon MPC to control the industrial distillation column. The controlled outputs are the level of liquid in the overhead drum (𝑦₁ %) and the temperature at tray number 68 (𝑦₂, °C). The manipulated variables are the steam flow rate to the reboiler (𝑢₁, t/h), the distillate flow rate (𝑢₂, m³/d), and feed temperature deviation (𝑢₃, °C). Figure 1 shows schematically the structure of the existing regulatory loops of the distillation column.

Figure 1 

Process diagram for alkylation system.

From practical tests, for a sampling period 𝑇=1 min, the following two models were obtained at different operating conditions: 𝐺₁(𝑧)=⎡⎢⎢⎢⎣2.3𝑧⁻⁴1−𝑧⁻¹−0.7×10⁻³𝑧⁻³1−𝑧⁻¹0.2𝑧⁻⁵1−𝑧⁻¹0.4604𝑧⁻⁸1−0.902𝑧⁻¹0.1915×10⁻³𝑧⁻³1−0.8632𝑧⁻¹0.03304𝑧⁻⁴1−0.9174𝑧⁻¹⎤⎥⎥⎥⎦(40)𝐺₂(𝑧)=⎡⎢⎢⎢⎣2.3𝑧⁻⁴1−𝑧⁻¹−0.7×10⁻³𝑧⁻³1−𝑧⁻¹0.2𝑧⁻⁵1−𝑧⁻¹0.1965𝑧⁻⁴1−0.9146𝑧⁻¹0.3215×10⁻³𝑧⁻²1−0.8852𝑧⁻¹0.03332𝑧⁻⁶1−0.9306𝑧⁻¹⎤⎥⎥⎥⎦(41) We observe that 𝑦₁ is integrating with respect to all the inputs while 𝑦₂ is stable. We should notice that there is uncertainty of about 20% in time constants and uncertainty of about 50% in process gains of the stable part of the model. Uncertainty in the dead times is also present in the model.

In the set point tracking problem simulated here, the desired value of the liquid level in the overhead drum (𝑦₁) was reduced by 1% and the desired value of the temperature at stage #68 (𝑦₂) was increased by 1°C. The tuning parameters of the controller that were used in the simulations included here are the following: 𝑚=2, 𝑄_𝑖=1, 𝑄_𝑠=1, 𝑅=diag([0.10.110]), 𝑆^𝑖_Δ=1×10³, 𝑆^𝑠_Δ=1×10³, and 𝑆^𝑖_𝑖=1×10³. Related to the values of the variables at the initial steady state, the input limits are 𝑢_max=[1040010], 𝑢_min=[−10−400−10], and Δ𝑢_max=[0.1500.01]. The nominal model is represented by model 𝐺₁(𝑠) defined in (40), and the true plant can be either 𝐺₁(𝑠) or 𝐺₂(𝑠) defined in (41). In this work, IHMPC represents the same controller as the one presented in Carrapiço et al. [13], which does not consider the uncertainty in the process model. In the first simulation, it is considered the case in which the true plant has the same model as the nominal plant. Figures 2 and 3 show the system responses for the nominally stable IHMPC proposed in Carrapiço et al. [13] and for the robust controller resulting from the solution to Problem 1 and designated here as IHRMPC. The two controllers have the same tuning parameters. One can see from Figures 2 and 3 that the nominal IHMPC and the robust IHRMPC perform similarly. This is quite surprising because it is widespread in the literature that the robust MPC should have a more conservative behavior than the controller based on the true model as this controller also takes into account the output predictions of model 𝐺₂(𝑠) that is quite different from the true model of the plant. This conservative behavior would produce a slower response, which is not observed in the simulation results obtained here. In the second simulation case, model 𝐺₂(𝑠) represents the true plant. The IHMPC is still based on model 𝐺₁(𝑠) while the proposed IHRMPC contains models 𝐺₁(𝑠) and 𝐺₂(𝑠). Figures 4 and 5 show the responses of the distillation system with each controller. We can see that the robust controller stabilizes the true model and the performance of the controller is acceptable. However, the IHMPC based only on the nominal model becomes unstable, which indicates that the controller based only on model 𝐺₁(𝑠) cannot be used to control the distillation column at the operating point where model 𝐺₂(𝑠) represents the true process.

(a)

(b)

Figure 2 

Outputs of the distillation system. Nominal case.

(a)

(b)

(c)

Figure 3 

Inputs of the distillation system. Nominal case.

(a)

(b)

Figure 4 

Outputs of the distillation system. Uncertain case.

(a)

(b)

(c)

Figure 5 

Inputs of the distillation system. Uncertain case.

5. Conclusion

In this paper, it was presented a new version of the robust infinitive horizon MPC with output feedback for systems with stable and integrating outputs. The adopted model formulation precludes the need to include a state observer, and the computer burden to run the controller may be reduced. To accommodate uncertainty in the process model, the state space model was built in two separate parts: one part is related to the integrating outputs and the other is related to the stable outputs. With this approach, it was possible to include the multiplant uncertainty in the model of stable outputs. The resulting optimization problem is a convex nonlinear program that can be easily solved with the available NLP packages. A simulation example shows that the implementation of the developed approach at real industrial systems may be achieved at least for systems of small to medium dimension.