Selection of Optimal, Controlled Variables for the ... - ACS Publications

Aug 5, 2010 - SINTEF, O.S. Bragstads Plass 2D, Trondheim N-7034, Norway. This article shows how controlled variables (CVs) of the regulatory control l...
0 downloads 0 Views 2MB Size
Download PDF
8624

Ind. Eng. Chem. Res. 2010, 49, 8624–8632

Selection of Optimal, Controlled Variables for the TEALARC LNG Process Finn Are Michelsen,* Berit Floor Lund, and Ivar Johan Halvorsen SINTEF, O.S. Bragstads Plass 2D, Trondheim N-7034, Norway

This article shows how controlled variables (CVs) of the regulatory control layer in a liquefied natural gas (LNG) plant can be chosen as linear combinations of measurements using self-optimizing control principles. By self-optimizing control, the CVs are chosen such that the set points of the CVs remain close to steadystate optimal despite disturbances, thus reducing the need for online reoptimization. Several methods for calculation of linear combinations within this framework are compared. Self-optimizing control design can also be used in the process design phase to place measurements by reducing a maximum candidate set of measurements to a best possible subset of measurements giving an acceptable loss. This article proposes a relatively simple method for successive selection (SS) of measurements and compares this approach to a more comprehensive branch-and-bound (BB) method for selection of measurements. The results indicate that, although the BB method gives lower average losses for very small subsets, the methods are comparable with respect to average losses for medium and large subsets, and the SS method outperforms the BB method in terms of computational load. 1. Introduction The choice of controlled variables (CVs) is an important step in control structure design,1 and it strongly affects the possibilities of obtaining optimal process operation. The best choice of CVs should be robust in the presence of unknown disturbances and process model uncertainties, and the final control implementation should be simple. In this article, the controlled variables of a liquefied natural gas (LNG) process are chosen as linear combinations of measurements using the self-optimizing control methodology.2-4 Self-optimizing control design is an approach to handle the important interaction between control structure design and process design. A set of controlled variables (CVs) is found that, when kept constant despite disturbances, uncertainties, and noise, provides acceptable losses in comparison with the use of an online optimizer that gives optimal operation. The set of CVs is found based on a profit criterion and constrained steady-state optimization, which normally is made as a part of process plant design. A simulation model of the TEALARC LNG process5 has been used as a case study for this work. This is a mixed-refrigerant process with two cooling cycles. Similar studies have been made of simpler refrigeration plants by Singh et al.6 and Jensen.7 Singh et al. studied an LNG plant designed and patented by SINTEF that is a single-refrigerant system with one compressor. In this article, several methods for calculating the set of CVs as proposed by Alstad et al.4 are compared, focusing on some of the most important disturbance variables. Also, it is shown how the results are affected by taking implementation error into consideration. In a process that is still in the design phase, the self-optimizing control methodology can be used to specify where to locate measurement sensors. Hence, an important part of the methodology is a search for a minimum set of measurements that should constitute the optimal CVs. Recently, Kariwala et al.8,9 suggested branch-and-bound methods for that purpose. Although other methods also exist for such search problems, this article compares one of their comprehensive methods with a relatively simple and computationally efficient successive selection method as suggested in this article. * To whom correspondence should be addressed. E-mail: [email protected].

In the following discussion, section 2 describes the objective function, disturbances, and constraints for the TEALARC LNG process. Section 3 describes the principles behind self-optimizing control design and the suggested search method in section 3.4.2. Section 4 reports results of the analysis. Section 5 analyzes and discusses the results. Conclusions and suggestions for further work are described in section 6. 2. TEALARC LNG Process and Model The purpose of an LNG plant is to turn natural gas into a liquid at atmospheric pressure for efficient storage and transportation by ships. This is done by first cooling the pressurized natural gas and then reducing its pressure. The resulting temperature for liquefied natural gas is normally about -162 °C, which is the condensation point of the lightest component, methane, at atmospheric pressure. In general, an LNG process setup consists of heat exchangers, compressors, separators, and valves. The TEALARC process was patented by TECHNIP10 and has been described by Paradowski and Dufresne and Paradowski et al.5,11 Based on those descriptions, the process used in this study is illustrated in Figure 1 and is the same process as studied by Wahl.12 The expander (Exp) and the storage tank (ST) are not included in the present model. Two refrigeration cycles, named the liquefaction and precooling cycles, are used in the plant. Both cycles contain refrigerants that are mixtures of components. The liquefaction cycle cools the natural gas in three heat-exchanger steps (HE1-HE3). This cycle contains two

Figure 1. Flow sheet of the TEALARC process used in this study.

10.1021/ie100081j  2010 American Chemical Society Published on Web 08/05/2010

Ind. Eng. Chem. Res., Vol. 49, No. 18, 2010

compressor stages, C1 and C2, with intercooling and aftercooling, as well as a separator. The mixed refrigerant used for the liquefaction of the natural gas is first partially condensed in heat exchangers HE4-HE6. This liquefaction refrigerant consists mainly of methane and ethane. In the separator, this refrigerant is split into a gas phase and a liquid phase with different boiling points. The liquid fraction, which has the highest boiling point, goes to the liquefaction heat exchanger (HE2). Thereafter, its pressure is reduced by the expansion valve V1. For an appropriate composition, this causes the temperature to decrease, which provides refrigeration to both the refrigerant and the natural gas. The gas fraction from the separator goes to the subcooling heat exchanger (HE3). The pressure and temperature are reduced by expansion valve V2 in order for the refrigerant to cool itself and subcool the natural gas. Precooling of the high-pressure liquefaction refrigerant is achieved by heat exchangers HE4-HE6. The precooling refrigerant is subcooled in HE4, where the pressure in a fraction of the stream is reduced in expansion valve V5. Further, the resulting partially vaporized stream is used to subcool itself and to precool the liquefaction refrigerant. Then, the vaporized refrigerant is recompressed. The remainder of the precooling refrigerant is then used in HE5. HE5 and HE6 are very similar to HE4, but the heat exchange takes place at lower temperatures and, thus, lower pressures for the precooling refrigerant. The recompressed streams from compressors C4 and C5 are externally cooled. 2.1. TEALARC LNG Model. There are some important model and simulator requirements for a self-optimizing control design analysis. First, a causal model is needed. This means that the model has to be formulated such that the effect of input changes, that is, manipulated and disturbance variables, on the measurements can be calculated. Second, the optimization calculations might have to handle discontinuities. Typically, commercial design simulators such as HYSYS are not able to handle these requirements.13 This issue is further discussed by Michelsen et al.,14 who developed a dynamic, control-relevant, causal model for the TEALARC LNG process. The model has a simplified two-phase flow description both in the NG and the refrigerants based on PT-flash calculations. Steady-state mass balances, dynamic energy balances, and average heat capacity values are assumed for the heat exchangers. The compressors are described by steady-state third-order compressor models. The control valves are described as expansion valves with estimated temperature and pressure drops. The model was implemented in Matlab and verified against a HYSYS model of the same plant, used for steady-state energy analysis. The nominal operating point generated by the energy analysis was also used for the dynamic model. Additionally, the parameters of the dynamic model were tuned to match the steady-state operating points to the HYSYS TEALARC model. 2.2. Operational Objective, Constraints, Degrees of Freedom, and Disturbances. Given a dynamic plant model dx ) f(x, u, d) dt y ) g(x, u, d, w) h(x, u, d) e 0

(1)

where x ∈ Rnx is the state vector, with nx as the number of states; u ∈ Rnu represents the manipulated variables (MVs), with nu as the number of MVs; d ∈ Rnd represents the unknown disturbance variables (DVs), with nd as the number of DVs; y ∈ Rny represents the measurements (Ys), with ny as the number of measurements; w ∈ Rny represents the measurement errors; f is

8625

the process model; h represents other constraints in the process variables including equality and inequality constraints; g is the measurement equation, that is, the relation between process states, disturbances, and sensor values. At steady state, dx/dt ) 0. Then fss(x, u, d) e 0

(2)

where fss includes f and h. The suggested function to be maximized in the operation of the process is the profit given by the difference between the income from the sale of LNG and the energy costs that are required to produce the LNG. This is the so-called objective function J max J (NOK/h)

(3)

J ) qlngVlng - qenergyVenergy

(4)

f(x, u, d) ) 0 y ) g(x, u, d, w) h(x, u, d) e 0

(5)

u

subject to

qlng is the production rate of LNG (kg/h); Vlng is the unit price of LNG (NOK/kg) (NOK denotes Norwegian krone); qenergy is the flow of consumed energy (kW); and Venergy is the unit price of energy (NOK/kWh). The objective function J is described in detail by Michelsen et al.14 Prices are assumed to be constant: Vlng ) 2.5 NOK/kg and Venergy ) 0.7 NOK/kWh. This means that the optimal operating point is given by maximum profit, which balances maximum income minus minimum energy costs. Important constraints in the process variables, h(x,u,d), are as follows: (i) 100% LNG fraction from the subcooler HE3, which can be controlled by the LNG temperature; (ii) dewpoint margin on the suction side of the compressors (e.g., 10 °C superheating),7 which provides safety constraints to avoid damage to the compressors; (iii) separator level (0-100%); (iv) valve opening (0-100%); (v) compressor capacity (maximum drive power and rotational speed, surge limitation); (vi) capacity in single components (e.g., maximum heat-exchanger duty); (vii) maximum feed (import limitation); (viii) maximum production (export limitation); (ix) storage capacity; (x) product specifications (composition and state); (xi) content limit of components that can freeze out in heat exchangers; and (xii) capacity limit of seawater cooling circuits. Degrees of freedom are the speeds of the five compressors C1-C5, the opening of the five expansion valves V1-V5, the flow rates of seawater to the four external coolers E1-E4, and the feed flow rate of NG. Jensen7 suggested that the amount of liquid in the cooling cycles, called the active charge, could also be considered as a manipulated variable for control. Figure 2 shows the profit as a function of the flow rate of natural gas, qNG, and the speed, NC2, of compressor C2. NC2 is given as a percentage of the maximum speed. All other degrees of freedom are constant. The optimum operation is located at 100% NC2. At this point, qNG is maximized, giving maximum production and maximum profit (i.e., at about 8800 kmol/h qNG). This is located at the steep edge of the profit surface. This edge constitutes a ridge of maximum points as a function of qNG at given NC2 values. For a given NC2, the profit increases gradually from low production rates until it drops steeply at a given rate because

8626

Ind. Eng. Chem. Res., Vol. 49, No. 18, 2010

Figure 2. Profit as a function of the flow rate of natural gas (qNG) and the compressor speed (NC2).

of refrigeration capacity limitations. Above this point, the refrigeration cycles are not able to reach down to the temperature of total liquefaction of the NG. For a given qNG value, the profit increases instantly from low compressor speeds in the range below liquefaction until a maximum point (threshold speed), where it decreases gradually for higher speeds because of increased energy costs for the cooling work. Hence, the optimum for qNG < 8800 kmol/h is not at maximum compressor speed. However, the degree of reduction in profit as a function of higher compressor speeds decreases with lower energy costs until no reduction at zero energy costs. In that case, the profit is given by the LNG income only, and the optimum is achieved at any compressor speed above the threshold speed. For combinations of low NC2 and high qNG, the refrigeration cycles are not able to reach down to the temperature of any liquefaction of the NG. In this range, the profit is negative and given by the energy costs.14 This means that the optimum operating point is given at maximum cooling capacity, which is given by maximum compressor speeds, valve openings, and cooling water flow rates in the external coolers, as energy consumption for the external coolers is ignored. This also means that the heat that can be removed from the NG stream is also maximal. Hence, cooling the LNG temperature below the dew-point temperature can only be achieved at the expense of reducing the throughput of NG. With these arguments, any of the degrees of freedom, except for the NG flow rate, used as a manipulated variable (MV) to control the LNG temperature at its dew-point temperature will be at its maximum constraint at the optimum operating point. Further, the optimum operating point at reduced production rate is also at the maximum constraint for the MV used to control the LNG temperature, which is obtained at reduced compressor speeds. This is shown in Figure 3 in the case of NC2 as the MV. NC1 is reduced to obtain optimum profit at 4 bar LNG pressure. The LNG temperature, dew-point margins, and separator level are commonly controlled to fulfill the constraints. Hence, in practice, seven of the degrees of freedom are used for these control loops. In this work, however, these variables are considered as uncontrolled constraints, which are verified to be fulfilled in the calculations based on constant values of all degrees of freedom that are not used for optimization. DVs can include exogenous changes affecting the system, process changes, changes in the specifications (constraints), and changes in the parameters (prices) that enter the objective function. They can also include parameter variations in the

Figure 3. Steady-state profit and control variables for controlling TLNG as a function of NC1.

process. Changes in heat transfer in the heat exchangers are regarded as the main disturbances. Such changes are typically caused by composition variations in the natural gas or in the refrigerants. Hence, these are the actual disturbance sources. The relation between such composition variations and the heat transfer in the heat exchangers is, however, not modeled. Instead, variations in the heat capacities of the natural gas and of the refrigerant in the liquefaction cycle are considered as unknown disturbances. Implementation errors can include control input deviations and measurement errors (MEs). Input deviations are not considered in this study. 3. Self-Optimizing Control Design The first local approach for selecting CVs within the selfoptimizing control framework was the approximate minimum singular value or the maximum gain rule described by Skogestad and Postlethwaite.15 Thereafter, several methods for calculation of the set of CVs were proposed in the literature.2-4,16,17 Most of these methods find CVs based on a common set of linear combinations of measurements and provide similar results.17 The methodology of self-optimizing control as explained in this article is based on the work by Halvorsen, Alstad, and Skogestad and co-workers.2-4 The basis for this methodology is to choose controlled variables such that changes in unknown (i.e., not measured or estimated) disturbances and other uncertainties give minimum operational losses. This is done through analysis of the objective function. At the nominal optimum, some process variables will normally lie at their constraints. A corresponding number of free variables are removed from the optimization problem, and the objective function is reformulated into an unconstrained optimization problem. The unconstrained, reduced objective function is then the basis for finding the controlled variables as a linear combination of measurements, c ) Hy. In the following sections, we show how the nominal steady-state optimum is found, how H is calculated in different ways depending on assumptions about the disturbance sensitivities and measurement errors, how the average loss is calculated, and finally how the search for a minimum set of measurements that constitute the optimal CVs is done. 3.1. Finding the Nominal Steady-State Optimum. The nominal steady-state optimum (u*,d*,y*) for the plant operation is found by minimizing a scalar objective function J with respect

Ind. Eng. Chem. Res., Vol. 49, No. 18, 2010

Figure 4. Feedback implementation of optimal operation. Process with disturbances d and measurement errors w.4

such that the map (transfer function) G ) HGy from u to c is square and such that control of c results in “acceptable” (i.e., optimal) operation of the process. The calculation of H depends on the number of Y values compared to the number of MVs and DVs, and the cases considered are discussed in the following sections. 3.2.1. Specific Number of Measurements ny ) nd + nu: Original Null Space Method. This is the case with zero disturbance sensitivity. Two methods for this calculation are as follows: (1) In this case, H can be obtained from a singular value decomposition of FT. Then, HF ) 0 or, equivalently, FTHT ) 0. Thus, selecting HT as the input singular vectors of FT, corresponding to zero singular values in FT, gives an orthogonal basis. F is the optimal sensitivity matrix of y with respect to the disturbances d, defined by

Figure 5. Combining measurements y to get controlled variables c.

to the MVs, u, which are available for optimization min J(x, u, d*) u

(6)

subject to eq 5. Before this optimization, the process is stabilized by using other available degrees of freedom that are not included in u. Further, no inequality constraints are assumed to be active at the optimum. Otherwise, a set of degrees of freedom (MVs) are separated from those that are available to fulfill the active constraints at the optimum. Moreover, the nominal measurement error is assumed to be equal to zero. Because the optimal u is a function of d, the unconstrained optimization problem can be written as min J(u, d) ) J[uopt(d), d] ) Jopt(d) u

(7)

For a disturbance d, eq 7 can be solved. This means finding the input vector u that minimizes the objective function for the particular d, as long as the set of active constraints has not changed. 3.2. Calculation of H. The linearized (at the optimum), local relation between u, d, and y is described by the matrices Gy and Gyd and deviation variables. A deviation from the optimum of the measurement vector y is then given by ∆y ) Gy∆u + Gyd∆d + w

c ) Hy

∂yopt i (d*) j F) ) ∂dj

[ ] dyopt 1 dd1

...

dyopt 1 ddnd

dyopt 2 dd1

...

dyopt 2 ddnd

l

dynopt y dd1

·

··

...

at d ) d*

(10)

l

dynopt y ddnd

where the superscript opt denotes the reoptimized measurements after perturbation by ddj of the disturbance dj. (2) An explicit expression for H is given by4 ˜ y-1 H ) Ho ) Mn-1J˜G

(11)

˜ y ) [Gy Gyd], and Mn } where J˜ ) [Juu1/2 Juu1/2Juu-1Jud], G -1 1/2 Juu (HGy) . Juu ) ∂2J/∂u2 and Jud ) ∂2J/(∂u ∂d) are the diagonal and offdiagonal submatrices, respectively, of the Hessian of J at the nominal optimal point (u*,d*). Note that the solution is not unique, as there are an infinite number of matrices H that satisfy HF ) 0. This stems from the freedom of selecting basis vectors for the null space. The linear model for the selected controlled variables can be written as

(8)

where ∆y ) y - y*, ∆u ) u - u*, and ∆d ) d - d*. The CVs, c ∈ Rnu (same dimension as the vector of unconstrained MVs), are represented by single Y values, or combinations of Y values, y. Linear combinations of the Y values are assumed

nu×ny

8627

(9)

. This set gives the smallest possible deviation where H ∈ R in the objective function when the disturbance changes. Thereby, their set points cs can be kept constant once optimized such that real-time optimization of these set points is not necessary; see Figure 4. The issue of selecting H can be viewed as a “squaring-down” problem, as illustrated in Figure 5. The number of output variables that can be independently controlled, nc, is equal to the number of independent MVs, nu. However, in most cases, the number of available measurements, ny, is larger than the number of independent MVs, that is, ny > nu. The issue is then to calculate the nonsquare matrix H, with dimension nu × ny,

∆c ) H∆y ) HGy∆u + HGyd∆d ) G∆u + Gd∆d

(12) and the degrees of freedom in the matrix H can be used to affect G ) HGy and Gd ) HGdy. Because H is nonunique, Mn can be chosen freely, for example, equal to the identity matrix I or such that G ) HGy ) I, which means that Mn ) Juu1/2. This latter choice gives a decoupled steady-state response from u to c. Note, however, that the product MnH is unique. The H matrices given by these methods have dimensions nu × nd + nu, such that HF ) 0, that is, H lies in the left null space of F. 3.2.2. Extra Measurements ny > nd + nu or Too Few Measurements ny < nd + nu: Extended Null Space Method. In many cases, the number of available measurements is larger than nd + nu. Then, information about measurement errors can be used to calculate c. The magnitudes of the disturbances d and measurement errors w are scaled by the diagonal matrices Wd and Ww, respectively

8628

Ind. Eng. Chem. Res., Vol. 49, No. 18, 2010

∆d ) Wddn w ) Wwwn

(13)

where dn and wn are assumed to be any vectors satisfying

|[ ]| dn wn

2

e1

(14)

Note that it is not possible to have zero loss with respect to implementation errors, because each new measurement adds a “disturbance” through its associated measurement error, w. The following two approaches to this issue are considered: (1) Optimal H for Combined Disturbances and Measurement Errors. In this case, the optimal measurement combination is obtained by solving an optimization problem formulated as H ) min ||M||F H

(15)

where || · ||F denotes the Frobenius (Euclidean) norm. M is given by M } -MnHF˜ ) -Juu1/2(HGy)-1HF˜

(16)

where F˜ } [FWd Ww] explains the term “combined” in this case. An explicit formula for H is given by using standard results for constrained quadratic optimization HT ) HcT ) (F˜F˜T)-1Gy[GyT(F˜F˜T)-1Gy]-1Juu1/2

(17)

This H matrix gives the average effect of combined bounded disturbances and measurement errors on the loss. Note that the matrix F˜F˜T needs to be full rank, which implies that this formula for H does not generally apply to the case with no measurement error, Ww ) 0, but otherwise, the expression for H applies generally for any number ny of measurements y. One special case, when the expression for H applies also for Ww ) 0, is when ny < nd, because then F˜F˜T remains full rank. This is also called the exact local method. (2) Optimal H for Disturbances Using Possible Extra Measurements to Minimize the Effect of Measurement Errors. An alternative approach is to first minimize the loss with respect to disturbances, and then, if there are remaining degrees of freedom, minimize the loss with respect to measurement errors. In this case, the optimal measurement combination is obtained by minimizing the error for the zero-disturbance sensitivity that leads to the explicit null space eq 11 ˜ y - J˜ E ) MnHG

(18)

An explicit formula for H in this case is given by ˜ y)†Ww-1 H ) He ) Mn-1J˜(Ww-1G

(19)

˜ y, and Mn are defined in the previous section. † where J˜, G denotes the pseudoinverse of the matrix. This H matrix minimizes the Frobenius (Euclidean) norm of E, ||E||F, and in addition minimizes the error sensitivity, ||Ww||F, among all solutions that minimize ||E||F.4 This minimization describes the average effect of the bounded disturbances and measurement errors on the loss. In the case ny < nu + nd, He is not affected by the measurement error. Then, eq 19 becomes ˜ †y H ) He ) Mn-1J˜G

(20)

Both of these approaches need the matrix Gy, which is calculated numerically by perturbing a dynamic model simulated to steady state. 3.2.3. Calculation of the Optimal Sensitivity Matrix F. F can be computed in two ways, which involve optimization for the selected disturbances using a nonlinear steady-state model of the plant: (1) numerical approximation of eq 10 by perturbing a dynamic model simulated to steady state and an optimizer and (2) calculation from the matrices Juu and Jud and steadystate gain matrices Gy and Gdy using F ) -(GyJuu-1Jud - Gyd)

(21)

3.3. Calculation of Average Loss. In self-optimizing control design, the loss is defined as the difference between J using constant set points and the optimal J. To calculate the former, a control structure with constant set points has to be implemented. Alternatively, the average loss can be calculated based on the above matrices. This was expressed by Kariwala et al.16 Lav )

1 ||M||F2 6(nu + ny)

(22)

Equations 17 and 19 are based on this approach. 3.4. Search Methods for Subsets of Measurements. If the self-optimizing control analysis is part of a design for a whole plant, an expansion of the method can be used to determine where to place measurements in the process. Given the ability with a process model to put sensors at any location in the plant or given the ability to put sensors at specific locations in the plant, it is valuable to identify the best available subsets of measurements that should be combined as controlled variables. Given Xm as an m-element set of all available measurements, the combinatorial optimization problem of subset selection involves finding an n-element subset, Xn ⊂ Xm, such that a selection (search) criterion is minimized among all possible Xn ⊂ Xm, as described by Kariwala and Cao.8 For small m and n, the globally optimal subset can be obtained through an exhaustive search. However, the brute-force method that solves this problem by enumeration of all m!/(m - n)!n! possible alternatives is usually not computationally useful for large-dimensional problems. Hence, effective search methods are needed that rapidly rank the importance of inclusion of the measurements. When plotting the loss against the number of included measurements, a gain of including extra measurements is most likely the largest when the measurement vector is small, but will decrease as measurements are included. This was also discussed by Alstad.18 Above a certain number of measurements, the gain in terms of reduced loss is most likely insignificant. Despite the fact that the analysis is local and based on linearization around a nominal point, the method can provide valuable input into determining where to locate actual measurements in the process. Such a method can be adjusted and extended by, for example, forcing the inclusion of some specific measurements if this is relevant. In this section, two such effective methods are presented: (1) partial bidirectional branch-and-bound (PBB) and (2) successive selection (SS). The principle differences between these methods are the search criterion and the search method. Specifically, for the search criterion, PBB finds the H matrix that minimizes the average loss (eq 22), whereas SS identifies the column of H with the least 2-norm. For the search method, PBB uses a partial bidirectional branch-and-bound technique, whereas SS employs iterative scanning of all columns of H. Both methods include

Ind. Eng. Chem. Res., Vol. 49, No. 18, 2010

all candidate measurements into a maximum measurement vector as the outset. 3.4.1. Partial Bidirectional Branch-and-Bound (PBB) Method. Kariwala et al.8,9 suggested branch-and-bound methods for the search for subsets of measurements whose combinations can be used as controlled variables. The partial bidirectional branch-and-bound method is especially applied to Hc (eq 17), which is valid for minimizing the average loss (eq 22). The methods are globally optimal, and they showed that the H matrix that minimizes the average loss is superoptimal in the sense that it also minimizes worst-case loss simultaneously.19 Cao20 published a Matlab code for this calculation. 3.4.2. Successive Selection (SS) Method. The principle of this method is that if one measurement enters each of the controlled variables with only a small factor, this measurement is assumed to be not important for the controlled variables. Hence, the corresponding row is removed from the matrices F, Gy, Gyd, and Ww. However, in doing so, H needs to be recalculated as the remaining elements of each row of H do not necessarily give minimum loss. The method can be summarized in this algorithm: 1. Start with a maximal measurement vector, with, for example, every element in the state vector of the model measured, and calculate the corresponding full matrices F, Gy, Gyd, and Ww. 2. While ny g nd + nu, do the following: a. Calculate H. b. Consider the columns of H. Find the column with the smallest 2-norm and remove the corresponding measurement from F, Gy, Gyd, and Ww, thus reducing the dimensions of these matrices (ny ) ny - 1). c. Calculate and save the average loss. 3. Repeat.

8629

Figure 6. Profit as a function of compressor speed.

4. Simulation Case 4.1. Case Description. In the analysis of the TEALARC LNG process, a large vector of candidate measurements is assumed, consisting of temperatures, flow rates, and pressures in every stream in the plant and the mole fraction of methane in the gas and liquid outlet streams from the separator. This gives a total of 66 possible measurements. As a result, the important variables to measure for control should be found, so that measurements with minor significance can be left out in the final implementation. The measurement vector is scaled by the assumed maximum values for the respective measurement value. An algorithm is implemented in Matlab for calculation of the matrices as described above. The DV and MV vectors d and u are perturbed by 0.1% in both the positive and negative directions from their nominal operating values. In all calculations, there is no liquid in the streams upstream of the compressors. As the focus of a self-optimizing control analysis is to choose controlled variables to counteract the effects of the disturbances, a key issue is, of course, which disturbances to include in the analysis. In the following sections, the results of a case with one MV and two DVs are shown. The calculations shown follow the steps in the procedure described in section 3. The speeds of the compressors are important MVs for LNG plants.7,21 Figure 6 shows the profit as a function of the speed NC2 of compressor C2 as the MV. This is a cross section of Figure 2. The optimum is located at a relative speed of about 68%, giving a maximum profit of about 328 kNOK/h. This is located at the steep edge of the profit curve. A consequence of

Figure 7. Average loss (kNOK/h) as a function of subset size of measurements.

the steep edge is that the Hessian is very large. This is in accordance with the results obtained by Singh et al.6 Unknown disturbances are described in section 2.2. There are no active constraints either at the nominal optimum or at the four optima related to perturbation of d (i.e., one for each of the two disturbances in the positive and negative directions). Small perturbations are chosen in order to make small linearization errors. The magnitudes of the disturbances d and measurement errors w as quantified by the diagonal scaling matrices in eq 13 are Wd ) 0.1%, and Ww includes 0.1% error for the temperature measurements (such sensors often have good accuracy) as well as the pressure, methane, and level measurements and 0.2% error for the flow rate measurements. 4.2. Results. 4.2.1. Comparison of Methods for Calculation of the H Matrix. To find a reasonable minimum number of measurements, Figure 7 shows the average loss, Lav, obtained using eq 22, as a function of the number of measurements based on Hc and He using the SS method described in section 3.4. Mn ) I for all calculations. The measurement selection is made by starting out with the vector of candidate measurements, that is, ny ) 66. Note that the equations are valid also for ny < nu + nd ) 3. However, the lower limit is ny ) nu ) 1. He is calculated according to eq 20 in these cases. Whereas eq 20 is valid only without measurement error, Hc is valid only including measurement error in these cases. Note also that Lav based on He (the

8630

Ind. Eng. Chem. Res., Vol. 49, No. 18, 2010

Table 1. Description of Measurements in the Optimal Measurement Combinations measurement number

description

1 5 12 40 41 51 56 57

natural gas temperature from cooler HE1 inlet refrigerant temperature of HE1 outlet temperature of E2 flow rate of refrigerant in cooler HE1 flow rate of liquid from the separator flow rate at the seawater side of E2 flow rate through compressor C3 flow rate through compressor C4

extended null space method) is larger than Lav based on Hc (the combined disturbances and measurements method). This result is related to the different ways in which measurement errors are handled by these methods. Figure 7 shows that Lav decreases considerably from ny ) 1 to ny ) nu + nd ) 3, that is, to the number of measurements corresponding to the zero-disturbancesensitivity case. A reasonable number of measurements in this case is five. By the same procedure, Ho is first calculated using the whole measurement vector and the pseudoinverse of Gy. Hence, by this method, the measurement vector has to be reduced to ny ) nu + nd ) 3 elements. Equations 23-25 show the optimal c calculated by the three different H vectors (H is reduced to a vector when nu ) 1) using the SS search method. cHo, cHc, and cHe are based on Ho, Hc and He, respectively. The H vectors are normalized by the Frobenius norm after reduction of the y vector.4

Figure 8. Average loss (kNOK/h) as a function of subset size of measurements.

cHc ) -0.38y12 - 0.35y40 - 0.35y51 - 0.59y56 - 0.51y57 (23) cHe ) 0.44y1 + 0.31y5 - 0.31y41 - 0.29y51 - 0.73y56 (24) cHo ) 0.62y1 - 0.74y41 - 0.26y51

(25)

cHc and cHe are different, although they contain some related variables. This means that the strategy of how to account for measurement error has some implications for the optimal measurement combination c. Further, all of the measurements related to cHo are included in cHe, as expected. This means that, because measurements 1, 41, and 51 constitute the optimal measurement combination c in the zero-disturbance-sensitivity case, measurements 5 and 56 are used to minimize the effect of measurement errors. The unscaled Ho and He matrices give HF ≈ 0; that is, H lies in the left null space of F. However, this is not the case for Hc because of the combined procedure of including the measurement errors. The sequential SS method gives a slightly different c than would be obtained by simply removing the 63 measurements in the original H matrix that have the smallest 2-norms. By following that procedure, measurements 40, 41, and 51 constitute the optimal measurement combination in the zero-disturbancesensitivity case. Table 1 describes these measurements. Note that the LNG temperature is not included in any of the optimal controlled variables. 4.2.2. Comparison of the SS and PBB Methods for Selection of Subsets of Measurements. Figure 8 shows the average loss, Lav, as a function of number of measurements based on Hc obtained using the SS and PBB methods described in section 3.4. The losses are about the same for the two methods for subsets larger than four. This confirms that a reasonable number of measurements in this case is five.

Figure 9. CPU time for the SS and PBB search methods.

Equation 26 shows the optimal c obtained from Hc and the PBB search method cHc ) -0.37y40 - 0.26y41 - 0.37y51 - 0.60y56 - 0.55y57 (26) Compared to eq 23, this c differs by only one measurement: y41 vs y12. This is also the case for larger subsets of measurements. Figure 9 shows the required CPU times for the two methods. For small subsets, the PBB method uses slightly less time than the SS method. For example, for a subset size of 5, the CPU times are 2.16 and 2.60 s, respectively, on a Windows XP SP2 notebook with an Intel Core Duo T7600 processor (2.33 GHz, 3.24 GB of RAM) using Matlab 7.1.0.246 (R14) SP3. However, for larger subsets, the SS method outperforms the PBB method. For the scanning of all subsets from 1 to 20, the CPU times are 322 and 6 s, respectively. Whereas the CPU time increases exponentially for the PBB method, it decreases about linearly for the SS method. 4.2.3. Verification of the Results Based on the Nonlinear Model. Because the self-optimizing control methods are local and linear, the loss is small only locally. Hence, it is essential to verify the results using the nonlinear model. Although the optimal c can be controlled by a proportionalintegral (PI) controller with zero steady-state control error at

Ind. Eng. Chem. Res., Vol. 49, No. 18, 2010

rejection of the heat capacity disturbances, there is a certain loss at the optimum operating point. By keeping c at its optimal set point, the loss is calculated as the reduction in profit from the optimum obtained at steady state by changing the two disturbances by (0.1%. The loss based on cHo is about 0.034% for the increase in the natural gas heat capacity and 0.032% for the decrease in the heat capacity of the refrigerant in the liquefaction cycle. The differences between the losses based on the full and reduced measurement vectors are negligible. Because liquefaction of natural gas is done by cooling the gas to the condensation temperature of the lightest component, which is methane, this LNG temperature could be expected to be an important variable to be controlled. Hence, for comparison with the self-optimizing control strategy, by control of the LNG temperature (i.e., a single measurement control strategy) the loss is about 0.035% for the increase in the natural gas heat capacity and 0.037% for the decrease in the heat capacity for the refrigerant in the liquefaction cycle, whereas the losses at no control are about 0.80% and 0.99% for these two perturbations. For the other two perturbations, the smallest losses are obtained at no control: 0.01% versus about 0.05% for the control cases. This is because the self-optimizing control methods are based on minimization of the average loss. Hence, better control strategies might exist for other disturbance cases. 5. Discussion The self-optimizing control design methodology includes model-based methods for control design. A process model, however, is used only in the design phase and not in the operation of the plant. Hence, gradual process model mismatch during operation, which does not affect structural properties that influence the optimal measurement combination, does not degrade the control performance. Further, the methodology is based on a steady-state, local, linear analysis. Hence, there is no guarantee that the controlled variables are globally optimal. However, they provide an indication of how sensitive the variables are to disturbances and measurement errors. The loss should be verified by simulations of the nonlinear model as in this study, and/or it should be better tested in the real (nonlinear) plant. Further, because of the steady-state assumption, if the set of active constraints changes, the set of optimal controlled variables has to be recalculated. This issue was considered by Cao.22 Alternatively, more advanced methods have to be applied. Caution concerning application of the results should, however, be used anyway because the thermodynamics is simplified. This leads to some additional uncertainty away from the nominal point. In practice, there is a need to include as few measurements as possible without too much loss in order to reduce the implementation costs. Hence, a key issue is reduction of the measurement vector with the objective of selecting a subset of the “best” measurements in the sense of those that have the largest influence on the optimal measurement combination. This can be done in several ways. In this study, the reduction is made by sequentially removing the measurement in the F matrix that has the smallest 2-norm in H and then recalculating H based on the new F matrix (SS) and by the branch-and-bound method by Kariwala and co-workers. By these methods, the most important measurements related to the influence of disturbance changes on the measurements are chosen. Alstad18 proposed an approach similar to the SS method based on the 2-norm of ˜ y ) [Gy Gyd], which chooses the most important measurements G related to the influence of nominal manipulated variables and disturbances on the measurements.

8631

Comparing the two search methods for subsets of measurements, the loss is slightly lower for the c vector based on the PBB method than for that based on the SS method, especially for small subsets. This indicates that the SS method is suboptimal, finding a slightly different subset of measurements for small subsets. This means that these methods are comparable with respect to reasonable sets of measurements being found for this problem. The difference between these methods regarding CPU time, however, varies significantly. As the SS method is based on successive removal of measurements from a candidate set, the CPU time decreases with decreasing subset size. For PBB method, however, the CPU time increases exponentially with increasing subset size. The methods were tested on a very small problem and found only one controlled variable. Because of the highly different relationships between CPU time and subset size for the two methods, a comparison for problems with larger dimensions should be carried out. The search methods are based on average losses. However, the results are qualitatively the same using the worst-case loss as described by Alstad et al.4 This was also found by Singh et al.6 for the SINTEF LNG plant. Proper scaling has to be made such that numerical errors do not corrupt the results. In this study, assumed maximum values for the measurements are used as scaling parameters. A reasonable alternative may be the normal variation range for each measurement. 6. Conclusions and Further Work This article shows how controlled variables (CVs) on the regulatory control layer in an LNG plant can be chosen as linear combinations of measurements using self-optimizing control principles. The results from a case study of the TEALARC LNG process indicate that the loss can be reduced by selecting the optimal set of controlled variables compared to a conventional approach using a constant-set-point control policy. This is shown for a large set and a reduced set of measurements including pressures, flows, and temperatures. This means that the optimal set of controlled variables can include several measurements; that is, controlling single measurements might not be optimal. Further, the optimal set of controlled variables is not necessarily located close to the manipulated variable for control, but generally has a low level of measurement error. Hence, there is a need for a systematic procedure for control design of this process, because it is difficult to obtain the best control structures for new plant configurations based on best practice and experience only. The objective function used to calculate the optimal solution for the case study of the TEALARC LNG process has a steep edge because of the cooling capacity limitations of the plant. This presents a significant challenge with respect to control. A small deviation in the production rate, caused by a process disturbance, will give a large reduction in profit if the disturbance is large enough. Hence, in this case, the operating point should be located slightly lower than the optimum, that is, at the lower-production-rate side. With respect to the selfoptimizing control design methodology, however, the calculations should be made at the optimum points. The original null space method is based on having as many measurements as the sum of the number of manipulated variables used for control and the number of disturbances considered for control design, that is, ny ) nd + nu. If extra measurements are available, that is, ny > nd + nu, then there are extra degrees of freedom in selecting H that should be used to reduce the sensitivity to measurement error. In that case, other

8632

Ind. Eng. Chem. Res., Vol. 49, No. 18, 2010

null space methods can be used. The combined measurement method appears to be the best of the two investigated methods in the sense of giving the lowest average loss. A suboptimal successive selection search method for finding a subset of measurements is comparable to a global optimal search method for the investigated problem. The methods investigated in this study find CVs based on a common set of linear combinations of measurements. As an alternative, Heldt17 proposed a method based on generalized singular value decomposition where individual measurement subsets are mapped onto each CV. This approach could be favorable in terms of flexibility, practical acceptance, and economic considerations. The self-optimizing control design methodology is not yet a mature technology. It is still under development, and there are several open issues to be studied. One issue is the robustness of the optimization calculation for larger numbers of variables, as this requires a large computational load. Another issue related to the optimization calculation is the convexity problem, which is always an important issue for nonlinear systems, although not in the case studied in this article. Further, this study has raised some new research questions. One question is determining which method for the selection of H is the best. Another is how to reduce the measurement vector with the objective of selecting the measurements that have the largest influence on the optimal measurement combination. Moreover, the implication of scaling of measurements should be investigated, as this should be done such that numerical errors do not seriously affect the results. A complete control structure design including all available manipulated variables for control also remains to be done. Acknowledgment This publication forms part of the Remote Gas project, performed under the strategic Norwegian Research program Petromaks. The authors acknowledge the partners Statoil, UOP, Bayerngas Norge, Aker Solutions, DNV, and the Research Council of Norway (168223/S30) for support. This article contains project information and preliminary results as a basis for final report(s). SINTEF accepts no responsibility of this article, and no part of it may be copied. Literature Cited (1) Skogestad, S. Control structure design for complete chemical plants. Comput. Chem. Eng. 2004, 28 (1-2), 219–234. (2) Halvorsen, I. J.; Skogestad, S.; Morud, J.; Alstad, V. Optimal selection of controlled variables. Ind. Eng. Chem. Res. 2003, 24 (14), 3273– 3284.

(3) Alstad, V.; Skogestad, S. Null Space Method for Selecting Optimal Measurement Combinations as Controlled Variables. Ind. Eng. Chem. Res. 2007, 46, 846–853. (4) Alstad, V.; Skogestad, S.; Hori, E. S. Optimal measurement combinations as controlled variables. J. Process Control 2009, 19, 138– 148. (5) Paradowski, H.; Dufresne, J.-P. Process analysis shows how to save energy. Hydrocarbon Process. 1983, 103–108. (6) Singh, A.; Hovd, M.; Kariwala, V. Control Variables Selection for Liquefied Natural Gas Plant. Presented at the 17th IFAC World Congress, Seoul, South Korea, Jul 6-11, 2008. (7) Jensen, J. B. Optimal operation of refrigeration cycles. Ph.D. Thesis, Department of Chemical Engineering, NTNU, Trondheim, Norway, 2008. (8) Kariwala, V.; Cao, Y. Bidirectional branch and bound for controlled variable selection. Part II: Exact local method for self-optimizing control. Comput. Chem. Eng. 2009, 33 (8), 1402–1412. (9) Kariwala, V.; Cao, Y. Bidirectional Branch and Bound for Controlled Variable Selection Part III: Local Average Loss Minimization. IEEE Trans. Ind. Inform. 2010, 6 (1), 54–61. (10) TEALARC. WIPO Registration Number 411865, TECHNIP, Courbevoie, France, 1974. (11) Paradowski, H.; Leroux, D.; Sguera, O.; Feresin, R. La Liquefaction des Gas Associes. Presented at the 7th International Conference on Liquefied Natural Gas, Jakarta, Indonesia, May 15-19, 1983. (12) Wahl, P. E. Enabling Production of Remote Gas: Pinocchio LNGsEnergy Analysis; Sintef Energy Research:, 2007. (13) Jøndal, H. Analyse av reguleringsstruktur for LNG-anlegg. Master’s Thesis, NTNU, Trondheim, Norway, 2008. (14) Michelsen, F. A.; Halvorsen, I. J.; Lund, B. F.; Wahl, P. E. Modeling and simulation for control of the TEALARC LNG process. Ind. Eng. Chem. Res., published online Jul 7, 2010, http://dx.doi.org/10.1021/ie901650e. (15) Skogestad, S.; Postlethwaite, I. MultiVariate Feedback Control. Analysis and Design; Wiley-Interscience: New York, 1996. (16) Kariwala, V.; Cao, Y.; Janardhanan, S. Optimal measurement combinations as controlled variables. Presented at the 8th IFAC Symposium on Dynamics and Control Process System (DYCOPS 8), Cancun, Mexico, Jun 4-6, 2007. (17) Heldt, S. On a new approach for self-optimizing control structure design. Presented at the 7th International Symposium on ADCHEM, Istanbul, Turkey, Jul 12-15, 2009. (18) Alstad, V. Studies on selection of controlled variables. Doctoral Thesis, NTNU, Trondheim, Norway, 2005. (19) Kariwala, V.; Cao, Y.; Janardhanan, S. Local self-optimizing control with average loss minimization. Ind. Eng. Chem. Res. 2008, 47 (4), 1150– 1158. (20) Cao, Y., Bidirectional Branch and Bound for Average Loss Minimization. Matlab Central, November 17, 2009; http://www.mathworks. com/matlabcentral/fileexchange/25870. (21) Singh, A.; Hovd, M. Dynamic Modeling and Control structure design for a Liquefied Natural Gas Process. Presented at the 2007 American Control Conference, New York, Jul 11-13, 2007. (22) Cao, Y. Direct and indirect gradient control for static optimization. Int. J. Autom. Comput. 2005, 2 (1), 60–66.

ReceiVed for reView January 13, 2010 ReVised manuscript receiVed June 17, 2010 Accepted July 2, 2010 IE100081J