An Optimal Feedback Control Strategy for Nonlinear, Distributed-Parameter Processes

: In this work, an optimal state feedback control strategy is proposed for non-linear, distributed-parameter processes. For different values of a given parameter susceptible to upsets, the strategy involves off-line computation of a repository of optimal open-loop states and gains needed for the feedback adjustment of control. A gain is determined by minimizing the perturbation of the objective functional about the new optimal state and control corresponding to a process upset. When an upset is encountered in a running process, the repository is utilized to obtain the control adjustment required to steer the process to the new optimal state. The strategy is successfully applied to a highly non-linear, gas-based heavy oil recovery process controlled by the gas temperature with the state depending non-linearly on time and two spatial directions inside a moving boundary, and subject to pressure upsets. The results demonstrate that when the process has a pressure upset, the proposed strategy is able to determine control adjustments with negligible time delays and to navigate the process to the new optimal state. and suggest


Introduction
In a bid to attain top efficiencies amid increasing competition, reduced profit margins, and rising concerns to protect energy and environment, modern day process operations mandate increasingly greater performance from control methods. This is a major challenge when dealing with the ubiquitous, non-linear and non-uniform processes, which are described by sophisticated mathematical models. In these cases, it is very difficult to obtain closed-form control laws and stability characteristics.
In general, the control of non-linear processes requires the solution of the associated optimal control problem [1][2][3]. This is done in non-linear model predictive control (NMPC) at each sample time of the process over a future time interval of fixed duration, or receding horizon. In a multi-layered scheme, real-time optimization (RTO) enables the determination of optimal set points to be tracked by faster acting closed loop control based on a simplified process model. These approaches provide a flexible structure to cover a wide range of optimality criteria [2][3][4]. In Refs. [5][6][7], it is shown that the knowledge of the necessary conditions of optimality (NCO) aids in devising efficient and robust real-time optimization schemes by utilizing the gradient information of the process under uncertainty.
However, in all aforementioned cases, there is always an inherent computational overhead associated with the optimal control determination, which is proportional to process non-linearity and non-uniformity. This overhead is often significant enough to cause feedback delays when the computation time exceeds the sample time interval. Such delays result in suboptimal solutions, or worse, infeasible operations [8]. Development of control techniques to surmount this challenge is an active area of research. Reviewed in 2016, the techniques include fast but typcially sub-optimal updates with negligible delays using a few iterations of the associated non-linear programs in serial, parallel

Theory
Consider the open-loop problem of finding the control u(t) at the optimum of the objective functional subject to the process model, or state equation where f depends on u (the control vector), y (the state vector), and its spatial derivatives; and p is a vector of functions representing boundary conditions. They depend on y at spatial boundaries, y ∂ , and its derivatives with respect to x (the vector of Cartesian co-ordinates x 0 , x 1 and x 2 ) and u. While y is constrained by Equation (2), u is undetermined over space and time. It may be noted that f represents a general, non-linear distributed-parameter process, and is not restricted to any particular structure. It is, however, assumed that the first and second order derivatives of f and F exist with respect to y and u. The problem is equivalent to finding the optimum of the augmented functional, where λ is the costate variable. Note that for simplicity, we have presented the objective functional in Lagrange form for fixed final time, and not included any (in)equality constraints since they could be handled by appropriately augmenting the associated Hamiltonian [16]. At the optimum of J, its variation, δJ = 0. Using the principles of variational calculus, this result eventually provides the necessary conditions-differential costate equation with zero-valued costates at the final time t f along with boundary conditions, and the stationarity condition, i.e., the variational derivative of J with respect to u being zero-which must be satisfied to obtain the optimal control policyû(t) and the corresponding optimal stateŷ(t).

Optimal Process Under Upset
For a process running with an optimal control policy, a process upset (i.e., a change in a parameter value) renders the state non-optimal. Given a sufficiently small upset, the change in J from the optimal state (ŷ,û) can be expressed as where δy = y −ŷ, and δu = u −û. Thus, the restoration of the process to optimality under the upset can be posed as the optimal control problem of finding the control adjustment δu that minimizes δ 2 J-the second variation of J-subject to the perturbed state governed by where ∂f/∂y is the total derivative taking into account the spatial derivatives of y as functions of y.

Incorporation of State Feedback
Discretizing the spatial domain using uniform grid point spacing ∆x i along each x i -direction, we obtain where is the Hamiltonian corresponding to the spatially discretized J in Equation (3) using L, M, and N grid points along, respectively, x 0 -, x 1 -, and x 2 -directions. While F ijk s are F values, λ and f are vectors of their respective elements at these grid points. Thus, for n state variables, Note that the minimization of δ 2 J, subject to Equation (4) with δu(t) as the perturbed control function, is the linear quadratic regulator problem of optimal control with finite time horizon [14]. In the necessary condition for the minimum of δ 2 J, we introduce the Riccatti transfromation, where δλ, y and δy have the same structure as that of λ in Equations (7) and (8), and is a block positive definite Riccati matrix, each element of which is a symmetric matrix of size LMN × LMN. Finally, we incorporate into the necessary condition, the feedback control law where K is the m × n gain matrix with elements similar to S i,j . With the above considerations, the necessary condition for the minimum of δ 2 J results in the differential Riccati equation where and a stationarity condition for δu that provides the equation for the gain, In the above equations, A and C have the same structure as that of S, while B and E have the same structure as that of K, and D is an m × m symmetric matrix.

Feedback Optimal Control Strategy
Based on the above development, following is the strategy to apply feedback control adjustment to a process susceptible to upsets in a parameter value:

Offline Calculations
For each, different value of a parameter expected to undergo upsets, determine and store a repository of the following: 1. the optimal stateŷ(t) by solving the open-loop optimal control problem, and 2. the gain K(t) by solving the differential Riccati equation and the gain equation

Real-Time Control
Execute the process with a nominal value of the parameter, as well as pre-determinedû(t). Whenever there is an upset in the parameter, do the following: 1. obtain from the repository theŷ(t) and K(t) interpolated at the current parameter value, and 2. obtain and apply the improved control which is the control law given by Equation (11).
It may be noted that the above strategy is different from model predictive control, which requires online solution of the optimal control problem at each sample time of a running process over a finite, receding time horizon. The past control approaches on distributed-parameter systems have been either open loop or focused on processes described by specific types of equations, such as linear, parabolic, hyperbolic, and elliptic partial differential equations [17][18][19][20][21][22][23].

Application
We apply the feedback control strategy proposed above to a distributed-parameter system of a labscale, cylindrical heavy oil reservoir subjected to gas injection for heavy oil recovery. The gas penetrates the vertical reservoir and reduces the viscosity of the oil, which in turn drains to the bottom and is recovered. The objective is to keep the oil recovery maximum despite disturbances in pressure (P) by adjusting the gas temperature versus time [T int (t)], which is the control function. Finally, we execute and validate the optimal feedback control strategy using computer simulation of the detailed process model, which was experimentally validated in an earlier work [24].
The upset considered for this process is pressure, P. Practical oil field applications often experience a phenomenon of forced fluid flow into the wellbore resulting from the deviation in mud hydrostatic pressure. This well control problem is termed as "kick" and it may result in severe operational difficulties [25]. A well kick needs to be properly compensated for so that it does not result in production irregularities.
Upsets in P affect the gas solubility (see Appendix A) at the cylindrical surface where T int (t) is controlled to help maintain oil recovery at the maximum level. Thus, the control function is a boundary condition for the equation of change of temperature in the system, which is simultaneously described by the equation of change of the gas mass fraction as well as the system height. While the system temperature and gas mass fraction vary with time and radial as well as axial directions, the height varies with time and the radial direction only.  Figure 1 shows the schematic of the experimental setup in which the gas (nitrogen) at a given pressure and temperature is injected to vessel containing the cylindrical heavy oil reservoir suspended from a load cell. With gas penetration, the viscosity of heavy oil drops, and it starts draining out from the bottom of the model into a collection funnel. The accumulated oil is directed through an online viscometer to a flash separation tank. The gas dissolved in the oil gets separated in a flash tank, and the gas volume is measured in two water columns. The gas-free oil from the flash tank is weighed and is the oil recovered from the reservoir.

Process Model
Based on first principles, the process model is given by the following three state equations: where T is the temperature inside the reservoir, w is the mass fraction of the dissolved gas in the reservoir, Z is the time-dependent reservoir height along the vertical direction, t is time, r is radial distance, and z is the vertical or axial distance. Moreover, α = K r Kρg where K r and K are, respectively, the relative permeability and permeability of the reservoir, g is gravity, ρ is oil density, and φ is reservoir porosity withĈ P as specific heat capacity and k as thermal conductivity.
It may be noted that D = D(w, P, T) is the dispersion coefficient of gas in heavy oil (refer to Appendix B), µ = µ(w, P, T) is the heavy oil viscosity (see Appendix C), andμ(r) is the average viscosity at z = 0, and is a function of r. Thus,μ(r) at any time is the average of the viscosity values at the corners [(r, 0), (r + dr, 0), (r, dz) and (r + dr, dz)] of a differential element at a given r and z = 0.
It is assumed that the flow of the live oil through the porous medium along the vertical direction is governed by Darcy's law, gas diffusion is dominant only along the radial direction, and the porous medium has constant α, ρ,Ĉ P , and k in the range of temperature variation [26].

Initial and Boundary Conditions
Initially, the reservoir is at room temperature T 0 with height Z 0 , and no gas mass fraction, except at the surface. Thus, initial conditions are: The boundary conditions are the specifications of T and w at the surface, and the radial symmetric conditions as follows: In the above equations, w int is the interfacial gas mass fraction (i.e., gas solubility) while T int is the interfacial temperature. The latter is a control function depending on time, and influences w int .
Equations (15)-(23) constitute the highly non-linear, distributed-parameter model of the heavy oil recovery process. This model is the basis for the optimal control problem that is formulated next.

Open Loop Optimal Control Problem
This problem involves the maximization of heavy oil production, i.e., the objective functional using interfacial temperature versus time, T int (t), as a control function. Here,ṁ(t) is the rate of recovery of the heavy oil from the bottom of the reservoir. The maximization of I is equivalent to the maximization of the augmented functional which is obtained by adjoining I to the state equations, Equations (15)- (17), with the help of the three undetermined costate variables-λ 1 (r, z, t), λ 2 (r, z, t), and λ 3 (r, t).

Control Function
It may be noted that the control function T int (t) is the temperature at: 1. the curved surface, i.e., T(R, z, t) for any height and time; 2. the bottom surface, i.e., T(r, 0, t) for any radial distance and time; and 3. the top surface, i.e., T[r, Z(t), t] for any radial distance and time.

Necessary Conditions for Optimality
At the optimum, δJ = 0. This gives rise to the necessary conditions that must be satisfied by the control function, T int (t), to be optimal. The conditions are the state equations, costate equations, and the stationarity conditions as follows.
Costate Equations These equations are: where for a property x,x denotes the average x in the differential element at a given r and z = 0; and subscripted r and z denote partial differentiation, e.g., T r = ∂T/∂r, T rr = ∂ 2 T/∂r 2 , etc.
The costate equations are subject to the following terminal and boundary conditions: λ 1 (r, z, t) = λ 2 (r, z, t) = 0 at Stationarity Condition This stems from zeroing out the variational derivative of J with respect to T int (t), and is as follows:

Computational Algorithm
The following computational algorithm was executed to determine optimal, open loop control policy at a given pressure: 1. Provide an initial guess for control values at each sample time instant in the interval. We applied the second order finite-difference approximations along the radial and axial directions to the partial differential, state and costate equations (PDEs). This step yielded a set of ordinary differential equations (ODEs) for each PDE at the grid points. The ODEs were then integrated using the fifth order Runge-Kutta-Fehlberg method with adaptive step-size control [27]. Grid-independence analysis was done to determine adequate number of grid points to be used in computations [24]. Broyden-Fletcher-Goldfarb-Shanno (BFGS) method [28] was utilized to improve the control values. The computational algorithm was programmed in C++, and executed on a desktop computer with Intel Core-I7 processor (64 bit, 3.60 GHz and 15.9 GB usable RAM). Appendix D provides the parameter values used in the computations.

Open Loop Results
The above algorithm was applied for five different pressures-24.5, 43.5, 50, 58.5, and 74.5 psi. In addition, operating temperature ranged from 25.0 • C to 90.0 • C. These operating ranges for pressure and temperature were selected based on the design limitations of the experimental apparatus used to validate the process model in previous studies [26,29,30]. At each pressure, the algorithm was executed several times with different initial guesses using 30 control values equispaced over the 163.9 min time interval. It was found that the algorithm converged to the same, final values of control versus time at a given pressure. Figure 2 shows the control values thus obtained. The corresponding values of optimal oil mass, and computation time and iterations taken by the algorithm are shown in Table 1.  It is observed that as pressure increases from 24.5 psi to 74.5 psi, value of optimal I increases from 19.3 g to 58.3 g, whereas the average optimal T int decreases from 77.9 • C to 56.6 • C. This result indicates that a higher pressure causes more gas to penetrate the reservoir and increase oil recovery. In addition, a lower gas temperature is more conducive to oil recovery at a higher pressure.
It may be noted that after quickly achieving a maximum, the optimal T int becomes mildly wavy in nature. Such behavior leads to temporary reversals in temperature and concentration gradients, which have been found to augment oil recovery in past studies [26,29,31].
Based on these optimal open loop results, we now implement closed loop control in presence of pressure upsets in the range, 24.5-74.5 psi.

Closed Loop Optimal Control
Following the development of Section 2, the Hamiltonian associated with Equation (25) is given by where the two states (T and w), costates (λ 1 and λ 2 ) as well as their spatial derivatives span the radial and axial directions, which are discretized using equispaced grid points indexed, respectively, by subscripts i and j. The third state (Z) and the costate (λ 3 ) span the radial direction only.

Control Methodology
For the online control of the process, we proceeded as follows. Right after the convergence of the open loop algorithm in Section 3.2.3 for a given pressure, the differential Riccati equation, Equation (12), was integrated based on the Hamiltonian given by Equation (34), and the gain matrix was computed from Equation (14). In this way, a repository of the optimal states and gains was obtained for five pressures, i.e., 24.5, 43.5, 50, 58.5, and 74.5 psi.
Next, the following online computational algorithm was executed for the closed loop control of the process, which was subjected to single or multiple step changes in pressure: 1. Begin the process simulation at a specified pressure, and apply the open loop, optimal control values determined earlier. 2. When there is or has been a pressure upset, replace the current control T int with where y is the current state while K andŷ are the gain and optimal state at the new pressure. The last two items are interpolated from the pre-determined repository of optimal states and gains.

Results
In this section, we present the results obtained from the application of the developed control methodology on the heavy oil recovery process subjected to one or more pressure upsets or step changes. Figure 4 shows the results for the closed loop control of the process after a step change of 17% in the initial pressure-from 50 psi to 58.5 psi-at 22.6 min. The open loop, optimal control policies at the five different pressures are also shown. It is observed that the control T int after the pressure upset approaches in about 27 min to the optimal policy corresponding to the new pressure. This adjustment leads to an improvement in the objective functional value, i.e., oil recovery. It gets raised from 42.6 g (at 50 psi without upset) to 51.0 g, which is close to 52.8 g (at 58.5 psi without upset). Figure 5 shows the results for a process run similar to the previous one but with the upset delayed to 73.5 min. It is observed that the control takes about the same time to approach the optimal policy for the new pressure, and yields 48 g of oil. This value is higher than that for 50 psi but lower than that in Figure 4. This is reasonable as the process is subjected to higher pressure (which facilitates oil recovery) but for a shorter duration though, due to the delayed pressure upset.   Figure 6 shows the results for the process initiated at 50 psi, and then subjected to five pressure upsets of 52.5, 56, 58.5, 62.5, and 67.5 psi at, respectively, 28.3, 39.6, 56.5, 73.5, and 84.8 min. It is observed that the control begins to decrease with the onset of pressure upsets, which are all positive. This is expected since lower control values are optimal at higher pressures. The oil recovery is 54.8 g, which lies between objective functional values of the neighborhood pressures, i.e., 52.8 g (for 58.5 psi) and 58.3 g (for 74.5 psi).

Tracking Property
To examine the tracking property of the control methodology, we consider the process control run of Figure 5, which is subjected to 17% pressure increment at 73 min. Figure 7 plots the reservoir temperature at the grid point: (i = 2, j = 3). It is observed that the temperature under feedback control progressively approaches the optimal reservoir temperature for that grid point at the new pressure.
For all grid points, Figure 8 plots the square of the error between the reservoir temperature under feedback control, and the corresponding optimal temperature at the new pressure. It is observed that the errors tend to insignificant values with time as the controlled state tracks the optimal state corresponding to the new pressure. The time derivative of the plotted error is of the order 10 −5 in the end. These results demonstrate that the developed control methodology has the capability to steer the process to the new optimum in case of a pressure upset.

Computational Overhead
With the modest computer configuration described in Section 3.2.3, it took less than a microsecond to carry out online computations to determine and apply a feedback correction to control.
Thus, the online feedback control has negligible time delay, which is very beneficial since it helps minimize the propagation of the effects of upsets on a process. As to storage requirements, it may be noted that the repository, which was determined offline, required a total space of just 0.185 MB.
Overall, the results indicate that the proposed control methodology has a potential for real-time, feedback control of non-linear and distributed-parameter processes. This methodology has the ability to take control action with minimal time delay, and track the optimal state. Of course, the field implementation would require suitable instrumentation including soft sensors for online state measurement/estimation at different grid point locations of the system.    Figure 8. Square of the error-between the feedback-controlled temperature, and the optimal temperature at the new pressure-at different grid points in the oil reservoir.

Conclusions
Based on the perturbation from an optimal state, a closed loop control strategy was developed for non-linear, distributed-parameter processes. The strategy incorporates offline computation of a repository of open-loop optimal states and gains corresponding to different values of a parameter susceptible to upsets. When an upset in the parameter occurs, control adjustments are quickly determined from the repository to attain the new optimal state. From the data of dissolved gas (nitrogen) in heavy oil, the gas solubility (w int ) was determined as a function of temperature and pressure. For this purpose, the selected operating pressures were 24.5, 43.5, 50, 58.5 and 74.5 psi. For each of these pressures, experiments were performed at 25, 50, 75 and 90 • C. Figure A1 shows the gas solubility thus obtained. It is observed that it rises with both temperature and pressure.

Appendix B. Dispersion Coefficient of Gas
The dispersion coefficient of gas (nitrogen) was determined as a function of its mass fraction in heavy oil utilizing a method developed elsewhere [30,32]. Shown in Figure A2, the dispersion coefficient increases with the mass fraction to a maximum and drops subsequently.

Appendix C. Heavy Oil Viscosity
The following correlation was used for the viscosity of heavy oil as a function of temperature ( • C), pressure (psi), and the gas (nitrogen) mass fraction: µ = a 0 + a 3 T + a 6 w P 2 + a 1 + a 4 T + a 7 w P + a 2 + a 5 T + a 8 w [g cm −1 min −1 ] The parameter values in the above equation are: