Model Predictive Control with Integrated Model Reduction for a Continuous Lactide Ring-Opening Polymerization Process.

Poly(lactic acid) production has received increasing attention, mainly due to its inherent biodegradable thermoplastic properties and to its renewable-resource-based composition. This process is affected by changes in the operating conditions and by raw material impurities which influence the reaction rate and degrade the polymer properties. As the system model is multivariable with coupled dynamics and constraints, linear model predictive control (LMPC) is employed here. A model reduction technique is proposed to obtain an approximate linear representation of the nonlinear system around the operating point to minimize the calculation cost of the controller. The proposed LMPC approach is validated by simulation and is compared to a proportional-integral controller and a nonlinear model predictive control. It is found that LMPC has a superior performance in terms of off-spec time when a disturbance occurs in the feed, and it can restore the target conditions better and faster.

Introduction

Whereas the vast majority of polymers are derived from either oil or natural gas, there is a growing trend toward producing them from renewable resources. In addition, many applications call for biodegradable and environmentally friendly polymers to minimize plastic waste, hence the importance of biodegradable polymers in the production of plastics.[1] Poly(lactic acid) (PLA) is a renewable thermoplastic made from renewable resources (sugar cane, corn starch, cassava roots, etc.). PLA is now a green alternative to nonbiodegradable polymers in a wide range of applications such as medical implants,[2] fibers,[3] and biomedical applications.[4] According to a recent review article,[5] PLA has multiple attractive properties such as very low cost, renewable character, high biocompatibility, excellent material properties, high transparency, and thermoplasticity, which open up a wide range of application fields. The two chemical methods used to produce PLA, melt polycondensation and ring-opening polymerization procedures, are discussed in this Article, along with the effect of several catalysts and the different polymerization conditions.

In the industry, the ring-opening polymerization (ROP) of the lactide monomer produces a high average molecular weight of the polymer, which is required in most applications.[6] However, this reaction is highly sensitive to impurities and variations in operating conditions. In particular, any impurities present will influence the final conversion and the average molecular weight of the polymer. Moreover, it is possible that the polylactide monomer itself contains unknown quantities of impurities.[7] Thus, it is crucial to use a control strategy for rejecting perturbations and ensuring the required quality and productivity of PLA production.

The classical PI controller was first used in 2016[8] to control a PLA process formed by two continuous reactors and a loop reactor. They considered the catalyst and cocatalyst flow rates at the inlet of the first reactor as control inputs, and they considered both monomer conversion as well as the pressure at the end of the single-loop reactor as outputs. In fact, the pressure has a direct correlation with the viscosity and thereby with the average molecular weight of the polymer. To deal with the high coupling between the inputs, the authors proposed an interesting way to modify the PI error (i.e., the difference exists between required and real outputs). A comparison of this strategy with the classical PI was performed and showed that the newly proposed strategy weakened the coupling and was able to overcome the in-feed disturbances.

However, It is important to note that the process is highly nonlinear and multivariable, with constraints on the inputs and outputs, and also, the control is only realized on the boundaries. A PI controller cannot handle all of these constraints, and therefore, a more sophisticated controller, such as a model predictive control (MPC), would be more appropriate and will thus be employed in this work to control the PLA process.

The MPC has been successfully operated in various fields.[9−11] It has proven its ability to cope with multi-input/multioutput (MIMO) processes and consider explicitly any possible constraints imposed on manipulated and/or tracked variables. Also, it is highly appreciated for time-delayed system control because of its predictive ability.

The MPC requires a reliable mathematical model that should be easy to manipulate and able to accurately predict the behavior of the plant. One of the most challenging issues of MPC is the resolution of the constrained optimization problem. For large dimensional nonlinear processes, the optimization problem, to be solved at every sampling time, may demand a very large computation time, without a guarantee to reach the optimum. Therefore, different linear approximations of the nonlinear model are often considered in order to formulate the overall task as a linear and convex problem. Thus, we speak about a linear model predictive control (LMPC) instead of a nonlinear model predictive control (NLMPC) for the case where a nonlinear process is used. For instance, linearized models were used within an MPC strategy to control a dentrification reactor,[12] building climate,[13] and a counter-current heat exchanger.[14] Laguerre functions can also be used to parametrize NLMPC, to reduce the simulation time and improve the optimization conditions.[15,16]

In this work, the multivariable control of a ring-opening continuous polymerization process of lactide is envisaged. This system is made up of three reactors (Figure ): two tubular ones on each side of a second loop reactor in the middle. The process model is constituted of sets of nonlinear partial differential equations (PDEs). The discretization in space of the PDEs representing the system results in a large number of nonlinear algebraic differential equations that increase the computational burden. Despite the good performance of the NLMPC that we used in our previous work,[17] it was much longer to run than the PI controller (depending on the disturbance level, the simulation time was between 760 and 2730 s in the NLMPC and about 138 s in the PI controller). Therefore, a reduced linear model is considered here and incorporated into the MPC strategy to minimize the calculation cost. The controller objectives are to maintain the product quality and the process productivity at the desired values, mainly in the presence of disturbances. The most important product property considered here is the average molecular weight of the polymer. That weight correlates with the medium loop’s viscosity and thus with the pressure drop. As control inputs, catalyst and cocatalyst flow rates are considered. Note that the control acts directly only at the entrance of the first reactor, while the outputs are measured just at the exit of the loop reactor (as all of the ingredients, including the control variables, are fed at the inlet of the first reactor). In fact, the most important conversion takes place in the loop reactor (70%), while the conversion in the last tubular reactor reaches 95%, which only slightly affects the polymer quality.

Figure 1

PLA process formed of two tubular reactors and a loop reactor.

The structure of this paper is as follows: In the Mathematical PLA Model section, we provide the nonlinear mathematical model of the PLA process. In the Controller Design section, we present a summary of the control objectives, along with a presentation of the principles of the three control strategies employed: (i) the novel LMPC strategy based on a discrete-time linear reduced state-space model, (ii) the NLMPC proposed in our previous work,[17] and (iii) the PI controller-based method proposed in ref (8). A comparison of the simulation results of the three used strategies is presented in the Results and Discussion section. The last section of the paper presents conclusions and perspectives.

Mathematical PLA Model

There are several sections and reactor types in the PLA process design (Figure ): a first prepolymerization tubular reactor (R1), a loop reactor (R2), and a tubular reactor (R3). By using the loop reactor (with a recirculation rate of r = 10, i.e., ratio between the recycled stream and the feed stream), the yield and the average molecular weight of the polymer can be increased. Both the mixing quality and the residence time in this reactor (“loop”) are similar to those of a continuous stirred tank reactor (CSTR), but with enhanced heat exchange.[18] Therefore, it will be modeled as a CSTR.

Material Balances and Average Molecular Weight of the Polymer

A model of an isothermal unit with flat radial concentration profiles is used to describe each tubular reactor, whereas the loop reactor is modeled as an isothermal CSTR (i.e., a single section). The reaction is performed in bulk (i.e., a homogeneous liquid phase composed of melted polymer mixed with liquid monomer). Based on the lactide ROP kinetic scheme,[8,19−21] the corresponding differential equations of material balances may be written as functions of the first four moments of the dormant (μi) and active (λi) chains as shown in Table (eqs 1–11). The reaction scheme and the development of the moments[22] are given in the Supporting Information in Tables S1 and S2, respectively.

Table 1

Corresponding Differential Equations of Material Balances of the PLA Process

From the moments of the active and dormant chain distributions, eq is used to calculate the average molecular weight of the polymer, Mw (“dead chains” have a negligible effect):MM represents the molecular weight of the monomer lactic molecule (half of the lactic molecule constitutes a monomer repeating unit).

Equation provides the adimensional pressure expression of the reactor, Φ:[23]L and X are the length of the loop reactor and the monomer conversion, respectively. η and ηref are the viscosity and the reference viscosity, respectively. p1 and p2 are parameters. Mw,ref, R, T, and Tref are the reference average molecular weight of the polymer, the universal gas constant, the temperature, and the reference temperature, respectively.

Initial and Boundary Conditions

Table summarizes the initial and the boundary conditions. Initially, the reactors are assumed to be full of monomer only. The boundary conditions include the inlet flow rates of monomer, catalyst, cocatalyst, and, as contaminants, acid and hydroxyl functionalities, introduced into the first reactor. We assume that both hydroxyl and acidic impurities can be present in the monomer, where the fraction of hydroxyl impurities to acidic impurities is equal to α (α = 0.5). Such impurities lead to the production of dormant species (eq 17).

Table 2

Initial and Boundary Conditions

Discretization of the Partial Different Equations

The discretization of the PDEs is done using the finite difference method. This method is widely used for the simulation of chemical processes due to its practical implementation.[24,25] The system of eqs 1–11 is implemented in Matlab that solves a finite difference scheme with Ng = 30 axial grid points for each tubular reactor.[8]

Controller Design

Control Objectives

The most important industrial requirements for PLA processes are predefined monomer conversion and average molecular weight of the polymer. Since the main production takes place in the loop reactor, both of these variables have to be regulated at the outlet of the loop reactor. The choice to use the outlet of the second loop reactor as controlled variables is due to the fact that the conversion is very low in the first reactor with low accuracy, while the outlet of the third reactor is far from the controller, with a very slow response.[8] As previously mentioned, these variables are strongly coupled since the increase in conversion correlates with an increase in the average molecular weight of the polymer (eq ). It is therefore important to consider this coupling in the controller design. Based on the kinetic scheme, such outputs can be effectively regulated through the use of both the catalyst (ppmcat) and the cocatalyst (meqROH) as manipulated variables (i.e., the inputs are U = [ppmcatmeqROH]).[8]

Regarding sensors, the on-line measurement of monomer conversion can be determined by near-infrared, Raman, or Fourier transform infrared spectroscopy, while the average molecular weight of the polymer is measured indirectly using the viscosity correlation (eq ). It is possible to determine the viscosity locally either via a viscometer[26] or by correlation with the pressure drop (eq ). As a result, the considered outputs of the PLA process are the conversion at the loop reactor outlet (y1 = X2) and the difference in pressure between the loop reactor inlet and outlet (y2 = ΔΦ2) (i.e., yM = [X2ΔΦ2]).

Model Predictive Control Structures: NLMPC and LMPC

Model predictive control (MPC) is an optimal control strategy that can be applied to both linear and nonlinear processes. This control approach requires a process model to predict variable evolution over a defined finite prediction horizon (NP). Thus, an optimal control policy involves a minimization of an objective function that includes generally one term related to the error between desired and real outputs and another term to penalize the effort of the applied control, all over a finite prediction horizon. The control action is optimized along a control horizon (Nc), which is shorter than or at most equal to the prediction horizon. Only the first determined control action is then applied to the system, and the optimization is carried out again at the next sampling time NS.

The MPC structure is shown in Figure ,[12] with the following components:

Figure 2

MPC structure.

the nonlinear process model given by eqs 1–13, including disturbances.

the prediction model (that can be a nonlinear model or a linearized process model, but without disturbances).

The MPC algorithm.

A, the reference model related to the two set-point signals Yset(k) = [X2(k) Φ2(k)]:

R, the reference model related to regulation:

Using the control structure shown in Figure , we can writeThe MPC controller’s objective is to keep the nonlinear process output yP at the set-point yA, and this allows us to bring the predictive model output yM to its set-point yd:Thus, the cost function, J(k), can be formulated as follows:[12,15,27]where Q is a positive definite matrix, and R is a positive semidefinite matrix. The second term of this cost function aims to limit excessive control actions.

A nonlinear prediction model of the process was employed in our work previously,[17] which leads to nonlinear MPC (NLMPC). At each sample time NS, the future sequence of manipulated variables was determined through solving the optimization eq (eq ). The main disadvantage of NLMPC is the complexity related with the resulting nonlinear optimization problem, which must be solved at every sampling time. Thus, it demands a very large computation time, and the convergence is not ensured for some model structures. Therefore, a reduced linear model is considered in this work and incorporated into the MPC strategy (LMPC case) in order to reduce the computational cost. In this case, an explicit solution of the problem can be found which reduces the calculation time.

Linear Model Predictive Control (LMPC)

Linear Model Reduction

A linear prediction model M is considered in this paper to describe the evolution of the nonlinear discrete-time state space model and thus to describe the influence of the two inputs variation (U = [Δu1Δu2] = [ΔppmcatΔmeqROH]) and the main disturbances (which are due to ΔmeqCOOH), on the two outputs variation ():where xm = [Δy1Δy2ΔmeqCOOH] ∈ R, U ∈ R, and yM ∈ R are the state variables, the inputs, and the outputs, respectively. AM, BM, and CM are state matrices with dimensions R, R, and R, respectively.

The nonlinear model of the process presented by eqs 1–11 is used to generate input–output data to allow an estimate of the model parameters of the two matrices AM and BM. This is done by varying the values of the inputs u1 and u2 and the main disturbances d around their nominal values.

LMPC Formulation

For the control strategy, the following notations are used (assuming Nc = NP = N):

y̆ = [y, ..., y], the outputs of the model A (eq ) along the horizon N

y̆ = [y, ..., y], the outputs of the model R (eq ) along the horizon N

y̆(k) = [y, ..., y], the outputs of the reduced model M (eq ) along the horizon N

ŭ(k) = [U(k), ..., U(k + N – 1)], the control signals along the horizon N

y̆ = y̆ – y̆, the desired outputs along the horizon N

The criterion in eq then becomeswhere , , and

yM(k) becomesThuswhereBy replacing y̆(k) with its value in eq , the following criterion can be obtained:[12]whereQ ≥ 0 and R > 0 imply φQ̆φ > 0 and . Thus, FU > 0, which guarantees the existence and uniqueness of the minimum of J(k). The solution of the problem that minimizes J(k) in the absence of constraintss on the input is given by

At each sampling time NS, the optimal inputs are calculated on the whole horizon, ŭ(k)*, but only the first element of this vector is applied to the process. Note that this calculation relies on constant matrices that only depend on the linear prediction model M. It is worth noting that we need a longer prediction horizon than the residence time of the PLA process (in the first tubular and loop reactor), i.e., NP > (tR2 + tR1). This is due to the distributed nature of the process, which means that this process responds with a delay to any input change. Thus, both LMPC and NLMPC will be inefficient if we vary the inputs without simulating the process for a significant time period, because the input variations will not have a relevant impact on the outputs. Therefore, for both LMPC and NLMPC, the prediction horizon NP is fixed at 40 min. Table summarizes the difference between the two strategies. The MATLAB optimization function “lsqnonlin” was employed in the NLMPC case to solve eq 31.

Table 3

Difference between the LMPC and the NLMPC

Control Strategy Based on a Modified PI

In this section, the strategy developed by Costa and Trommsdorff[8] is recalled to underline the shortcomings of its use and therefore to emphasize the necessity of the LMPC approach. The principle is to use the classical PI controller with a modified formulation of the error to account for inputs–outputs coupling (Figure ).

Figure 3

Modified PI version used in the work of Costa and Trommsdorff.[8]

The following equations show the differences that were employed in the modified PI version:

By looking at eq , ΔΦset is the desired set-point only if X2 equals X2set; otherwise, the current set-point is dynamically adjusted whenever X2 is different from X2set. This section presented the principles of the three control strategies employed to control the PLA process. The PI strategy is easy to develop and fast to compute, and an interesting modification was suggested to weaken the coupling between the inputs and outputs;[8] however, it cannot handle all of the process nonlinearity and constraints. Model predictive control can handle coupling between the inputs and outputs, process nonlinearity, and constraints, but it requires a precise process model. First, we proposed to employ NLMPC, where the criterion is straightforward to develop and can be computed using available optimization routines. We expect the computation time to be longer than the PI controller, so care should be taken to ensure that it is possible to realize on-line. Indeed, the nonlinear optimization task is usually nonconvex, for which the optimal solution is not guaranteed, and it may cause a high computational cost. To cope with this disadvantage of NLMPC and improve the optimization performance, LMPC was proposed. The LMPC formulates the overall optimization task as a linear and convex problem. However, this requires a further calculation to linearize the model around the interesting operating points. Also, note that it is inappropriate to use linear models to represent a very highly nonlinear model.

Results and Discussion

Simulation Conditions

The parameter values of the model are given in Table .[8,17,18,21] The set-points are set as the desired monomer conversion, X2set, at the loop reactor outlet equal to 0.6887 and the final desired average molecular weight of the polymer, Mw3,out, at the third continuous reactor outlet equal to 2 × 105 g mol–1. These two conditions result in a pressure drop, ΔΦset, in the loop reactor equal to 0.1419.

Table 4

Parameters Values of the ROP Process

In the following section, we will plot the figures as a function of τ, which is the adimensional ratio between the actual time t and the residence time in the loop reactor, i.e., τ = t/tR2. For all of the scenarios, at first, we simulate the PLA model over 20τ until the nominal steady state is reached, and then, a step-change in the acidic impurity level ΔmeqCOOH is made in the feed, which mimics a sudden impurity perturbation in the monomer concentration (positive or negative).

The LMPC controller is compared to the modified PI controller (eqs and 33)[8] and to the NLMPC (eq 31).[17] In the case of the modified PI controller, we assume a perfectly measured disturbance and that the computed inputs are implemented each minute (Ns = 1 min). However, for the cases of the LMPC and the NLMPC, both of these methods are not based on any measured disturbance. Thus, first, a method is developed in order to estimate the disturbances which can be then implemented for improving the LMPC and the NLMPC performances. Both the LMPC and the NLMPC suppose that the outputs are measured on-line and apply a new control value each 5 min (Ns = 5 min), which is more realistic.

Estimation of the Disturbance

For both LMPC and NLMPC cases, the level of the disturbance is first estimated and implemented in each strategy. Therefore, the nonlinear process model is employed to generate input–output data while changing the disturbance values ΔmeqCOOH as input in positive and negative directions and by recording the outputs variation of X2 and ΔΦ2. The ΔΦ2 variation was found to be more significant than the X2 variation when ΔmeqCOOH was varied and was therefore adopted for the estimation of the perturbation. The MATLAB routine “cftool” was used to identify the two following equations which give the relation between ΔmeqCOOH and ΔΦ2 in both positive (eq ) and negative (eq ) cases:These two estimated equations were then used for estimating the disturbance once a change in the outputs occurs, which then was implemented in the LMPC and the NLMPC strategies. It is assumed that, only after 15 min of its appearance, the disturbance can be estimated for both LMPC and NLMPC strategies. This is because, before 15 min, no impact on the output slopes could be seen. Thus, during the first 15 min, u1 and u2 were fixed at their nominal values. Then, from t = 15 min, the NLMPC optimization (eq 31) and the LMPC based on the identified linear model are performed using the estimated perturbations.

Identification of the Reduced Model Parameters

In order to identify the parameters of the linear model (eq ), the nonlinear process model represented by (eqs 1–11) was simulated using the nominal feed values over 20τ. Then, multiple stepwise changes were realized in the two manipulated variables and the acidic impurity level, around the nominal values to better capture the behavior of the two outputs.

The generated input–output data were used to estimate the two state matrices AM and BM of the linear model. This was done by optimizing the error between the outputs obtained from the nonlinear and linear models, using the routine “fminunc” of MATLAB. The obtained AM and BM are as follows:

Comparison between the PI, the LMPC, and the NLMPC

Various levels of positive and negative perturbations were supposed to occur in the inlet flows.[17] We added positive perturbations in the range 1–10 mmol kg–1 together with different negative perturbation levels (−1 to −5 mmol kg–1) to the nominal value which is equal to 5 mmol kg–1. The responses of the PI controller, the NLMPC, and the LMPC strategy to these perturbations were studied.

The PI tuning parameters were determined following the Cohen–Coon tuning method that was found to give a better performance compared to other methods, like Tyreus–Luyben and Ziegler–Nichols.[8]

Response to Positive Disturbances

Various levels of positive perturbations are supposed to take place in the input flows. Figure shows a comparison of the monomer conversion results, X2, provided by the modified PI, the LMPC, and the NLMPC. From these simulations, both the LMPC and NLMPC strategies appear to perform better than the PI controller. They are able to reject disturbances better and reach the nominal state faster, even without a direct measurement of the disturbances. Both strategies LMPC and NLMPC manage to reach the set-point X2 at the same time. The output X2 obtained by the LMPC strategy has a slightly lower peak than the one produced by the NLMPC strategy X2.

Figure 4

Time evolution of the X2 for different levels of positive disturbances.

Regarding the results of ΔΦ2 obtained by the PI controller, the NLMPC and the LMPC are compared in Figure for the same positive impurity levels. This figure demonstrates that the higher the disturbance, the bigger the pressure drop impact, ΔΦ2, as well as the more time it takes to reach the nominal state. The LMPC strategy performs much better than the PI controller and the NLMPC strategy in this case, as it converges faster than they do with a lower deviation.

Figure 5

Time evolution of the pressure drops ΔΦ2 for different levels of positive disturbances.

Response to Negative Disturbances

The response of the PI controller, the NLMPC, and the LMPC to a variety of negative perturbation levels was also considered. In this case, Figure shows that the three proposed strategies have an almost similar behavior in terms of X2 monomer conversion. All strategies can reject the perturbation within a reasonably short time.

Figure 6

Time evolution of X2 for different levels of negative disturbances.

The results of ΔΦ2 obtained by the three strategies for the same negative impurity levels are compared in Figure . The drift of the modified PI controller is clearly larger than that of the NLMPC and LMPC strategies. The LMPC strategy again gives better behavior than the PI controller and the NLMPC, with a lower peak than the other two strategies.

Figure 7

Time evolution of ΔΦ2 for different levels of negative disturbances.

Another advantage of using the LMPC is that the time to execute is low compared to the NLMPC strategy [LMPC (between 27 and 32 s), NLMPC (between 760 and 2730 s)]. It appears that the NLMPC takes more time to execute than the LMPC and the PI controller (between 138 and 140 s); this is due to the optimization routine which is carried out at each sampling time. Thus, there is the need to apply the LMPC to reduce the simulation time, and it is therefore possible to perform it online.

Conclusions

In this work, three control strategies were employed to control a highly nonlinear PLA process. This reaction is sensitive to impurities (mostly present in the monomer feed), so it is necessary to develop a control strategy for recovering the nominal operating conditions when any disturbances occur. The objectives of the control are linked to the productivity of the process (monomer conversion) and to the quality of the product (average molecular weight of the polymer). These outputs may be controlled through the manipulation of both catalyst and cocatalyst feed rates. Since the inputs/outputs are highly coupled, and the system is subjected to different constraints, a more advanced controller like NLMPC and LMPC is required. These two strategies were compared to a PI controller which was modified to deal with the coupling.[8] The LMPC was found to recover nominal steady states earlier than both the PI controller and the NLMPC. The use of a reduced model instead of the complex model, which is the case for the NLMPC, allows a reduction in the simulation time and thus the possibility of the implementation of the LMPC on-line. In this process, a disturbance estimate is needed to improve the LMPC behavior. Since many types of disturbances can occur, it is difficult to measure them. Thus, a simple estimated equation was used to approximate these disturbances. Even if this estimated disturbance is not exactly the same as the real one, the LMPC manages to restore the nominal conditions.