IGCC Power Plant Dispatch Using Infinite-Horizon Economic Model

Article pubs.acs.org/IECR

IGCC Power Plant Dispatch Using Infinite-Horizon Economic Model Predictive Control Benjamin P. Omell and Donald J. Chmielewski* Department of Chemical and Biological Engineering, Illinois Institute of Technology, Chicago, Illinois 60616, United States ABSTRACT: The Integrated Gasification Combined Cycle (IGCC) possesses many benefits over traditional power generation plants, ranging from increased efficiency to flex-fuel and carbon capture opportunities. A lesser-known benefit of the IGCC configuration is the ability to load track electricity market demands. The idea being that process modifications to enable dispatch capabilities will allow for a time-shift of power production away from periods of low energy value to periods of high value. The work begins with an illustration of Economic Model Predictive Control (EMPC) as a vehicle to exploit dispatch capabilities by pursuing the objective of maximizing revenue directly. However, implementation of EMPC can result in unexpected and, at times, pathological closed-loop behavior, including inventory creep and bang−bang actuation. To address these issues, an infinite-horizon version of EMPC is developed and shown to avoid many of the performance issues observed in the finite-horizon version. The paper concludes with an in-depth discussion of energy value forecasting and how the quality of forecasts can be incorporated into the design of the infinite-horizon EMPC controller.

I. INTRODUCTION Power plants based on the Integrated Gasification Combined Cycle (IGCC) possess many benefits over traditional power generation, including increased efficiency and flex-fuel and carbon capture opportunities. In short, a conventional IGCC power plant operates as in Figure 1, but without the storage units or MeOH plant. The gasification block converts coal and oxygen into synthesis gas, which is subsequently cleaned and decarbonized to output a stream of nearly pure hydrogen. This hydrogen is then converted to electric power by a combined cycle power block. A portion of generated power must be used to drive compressors for air separation and CO2 sequestration operations. Another possible benefit of IGCC is smart grid coordination through power output dispatch. Specifically, it is postulated that an IGCC plant can change power output to track grid demand and thus exploit electricity market price fluctuations. The ability to change power output from an IGCC plant can be achieved by a variety of hardware configurations (see Figure 1). The configuration most commonly cited is poly generation. In this case, a portion of the hydrogen from the acid gas removal unit is diverted elsewhere during periods of low energy value. This diverted hydrogen could be sold directly or chemically converted to a liquid fuel (possibly MeOH) for transportation applications.1 Another possibility is to store the hydrogen during periods of low energy value, and then draw from the inventory during peak demand periods. A similar scenario, with respect to the Air Separation Unit (ASU), is also possible. In this case, a compressed air storage unit is placed between the air compressor and the cryogenic distillation units. Such a configuration would allow for a constant flow of compressed air to the cryogenic distillation, while allowing the power to the air compressors to be anticorrelated with energy value, which would add power to net output during periods of high value. Similarly, significant power is consumed to compress CO2 to sequestration pressure. Thus, part of the captured CO2 can be stored at an intermediate pressure during periods of high energy value and eventually © 2013 American Chemical Society

pressurized during low-value periods. It is also noted that an unmodified plant has dispatch capabilities, in the sense that the gasification block can be asked to change hydrogen production rates. However, the ramp rate of the gasification block is expected to be quite slow, mostly because of the slow response time of the cryogenic distillation portion of the ASU.2 A fundamental question associated with smart grid coordination is the development of a control structure capable of accepting information from the smart grid. The first issue concerns set-point tracking abilities. Specifically, one must design regulatory level loops to achieve desired set-point ranges and ramp rate capabilities. Several efforts have explored this issue for a variety of IGCC configurations, as well as specific unit operations.1−4 While controller design at the regulatory level is clearly an important activity, such efforts typically lack sufficient motivation. That is, the performance objectives of desired setpoint ranges and achievable ramp rates are set somewhat arbitrarily. To arrive at sufficient motivation for these objectives, one must consider the supervisory controller, which will ask the regulatory controllers to implement its set-point commands. Two scenarios are possible. The first is a situation in which the regulatory loops are overdesigned. That is, the controllers are capable of achieving performance greater than the supervisory controller would ever request. In this case, there is likely overexpenditure regarding control system hardware. The second scenario is a situation in which the regulatory loops are incapable of achieving the performance desired by the supervisor. In this case, the regulatory loops could be redesigned (by adding more effective hardware) or the performance limitations of the regulatory loops could be communicated to the supervisor, so Special Issue: Process Engineering of Energy Systems Received: Revised: Accepted: Published: 3151

April 5, 2012 January 7, 2013 January 9, 2013 January 9, 2013 dx.doi.org/10.1021/ie3008665 | Ind. Eng. Chem. Res. 2013, 52, 3151−3164

Industrial & Engineering Chemistry Research

Article

Figure 1. Simplified process diagram of an IGCC power plant (with carbon capture) and possible equipment upgrades to achieve dispatch capabilities.

II. ECONOMIC MPC Consider a process model: ṡ = f(s,m,p), q = h(s,m,p), where s is the state vector, m is the vector of manipulated variables, p is the disturbance vector, and q contains the performance outputs. In addition, assume the following restrictions on the performance outputs: qmin ≤ q ≤ qmax. To implement a predictive-type controller, first, the continuous-time model must be converted to discrete-time form and, second, the notion of a predictive time index must be introduced. The resulting model is

that it knows the level of performance that is possible and can plan accordingly. One way to define a supervisory controller is to assume a merchant style of operation. Specifically, set the goal to be maximization of average revenue, defined as the integral of the product of energy value and net power produced. In this case, the energy value signal would be an input to the controller. However, unlike a disturbance, the objective is not to kill its impact on the output. Similarly, unlike a set-point command, the objective is not to force the output to track energy value (however, as we will see in subsequent sections, the set-point analogy is very close to the resulting objective). One way to approach this problem is through the use of Economic Model Predictive Control (EMPC).5 In this case, the quadratic objective function typically used in MPC6 is replaced by an expression directly reflecting revenue. Similar approaches have worked well in the area of heating ventilation and air conditioning (HVAC) control,7−12 as well as in the area of process operations scheduling.13−16 In both of these applications, computational issues, resulting from the large size of the online optimization problem, have been reported. Efforts to reduce this computational burden, by reducing prediction horizon, have resulted in undesirable closedloop characteristics, including inventory creep15 and bang−bangtype actuation.12 It is also noted that efforts to ensure closed-loop stability of the EMPC algorithm have required the use of specialized analysis techniques and nonintuitive additions to the basic algorithm.17−19 The next section outlines the EMPC approach and illustrates, through a simplified IGCC example, its ability to maximize revenue. In addition, the closed-loop eccentricities of the EMPC are illustrated. Section III introduces a new class of linear controllers. This new class, denoted Economic Linear Optimal Control (ELOC), is shown to possess a policy similar to EMPC, but only capable of enforcing equipment constraints, from a statistical perspective. In Section IV, the infinite-horizon EMPC controller is developed using ELOC as a basis. Section V discusses the important issue of forecasting energy value and the development of a measurement based formulation of the infinite-horizon EMPC policy.

s(k + 1|i) = fd (s(k|i), m(k|i), p ̂(k|i))

(1)

q(k|i) = hd(s(k|i), m(k|i), p ̂(k|i))

(2)

q min ≤ q(k|i) ≤ q max

(3)

s(i|i) = s (̂ i)

(4)

where the index i represents the actual time of the process and the index k is the predictive time. The idea being that at a time i, the controller will be given an estimate of the initial condition, ŝ(i), along with forecasts of the process disturbances, p̂(k|i). The controller then must select a sequence of control inputs m(k|i), k = i, ..., i + N − 1, such that the inequalities of expression (3) are satisfied. Since there may be many sequences satisfying expression (3), an additional performance measure is applied via the following optimization problem: ⎧i+N−1 min ⎨ ∑ g(s(k|i), m(k|i), p̂(k|i)) s(k | i), m(k | i)⎪ ⎩ k=i ⎪

⎫ ⎬ + gN (s(i + N |i), m(i + N |i), p ̂(i + N |i))⎪ ⎭ ⎪

(5)

subject to s(k + 1|i) = fd (s(k|i), m(k|i), p ̂(k|i)) q(k|i) = hd(s(k|i), m(k|i), p ̂(k|i)) q min ≤ q(k|i) ≤ q max s(i|i) ≤ s (̂ i) 3152

dx.doi.org/10.1021/ie3008665 | Ind. Eng. Chem. Res. 2013, 52, 3151−3164


Article

Figure 2. Simplified process diagram of IGCC power plant with H2 storage.

If the solution to expression (5) is m*(k|i), k = i, ..., i + N − 1, then the input given to the process s(i + 1) = fd(s(i),m(i),p(i)) is m(i) = m*(i|i). Then, at the next time step, i + 1, new ŝ(i + 1) and p̂(k|i + 1) estimates are calculated (based on new measurements) and given back to expression (5) to calculate m(i + 1) = m*(i + 1|i + 1). In the traditional implementation of MPC, g(·) and gN(·) are selected to be a quadratic functions with minimum values occurring at the target operating condition, where the target operating condition is selected prior to the MPC implementation. These smooth convex functions provide computational advantages, as well as a good starting point for stability analysis. In the case of EMPC, g(·) is selected to represent operating cost of the process. In many cases, the appropriate function is a simple linear relation. Since a linear function has no minimum, selection of the constraints, qmin ≤ q(k|i) ≤ qmax, is critical to the definition of a meaningful controller. Example 1: The operating conditions of Case 2 of the NETL Baseline Report20 form the basis of all examples. Specifically, the following conditions were assumed for the no dispatch design: nom nom nom Pnom N = 544 MW, PG = 642 MW, PAC = 67 MW, PCC = 31 MW, nom nom nom νH2 = νG = 90 tonne H2/h, νC = 790 tonne compressed air/h and νnom CO2 = 457 tonne CO2/h. For the current example, assume only two upgrades are performed (see Figure 2). The first is to add a hydrogen storage unit, capable of holding 700 tonne of H2. The second is to double the size of the power block, so that its maximum output is 1284 MW. Furthermore, assume that only these two units are capable of dynamic operation and all other units remain at the nominal operating conditions. Remark: Note that the equipment upgrades of this example are selected arbitrarily and intended only to illustrate the EMPC procedure. However, it is expected that changes to the assumed equipment upgrade sizes will impact the revenue increase expected from EMPC operation. Specifically, larger upgrades are expected to yield larger revenue increases. It is additionally recognized that a proper economic assessment of the IGCC dispatch opportunity will require a net present value calculation that weighs revenue gain against the capital cost of equipment upgrades. While the preliminary results of Yang et al.21 propose such a procedure, many of the engineering challenges (associated with H2 storage facilities and equipment degradation resulting from cycling of the power block) have yet to be addressed.

In sum, the current effort has the limited scope of assessing revenue gains under the assumption of prespecified and technologically viable equipment upgrades. As a final point, the assumption of lower bounds at zero for both storage capacity and power block throughput is simply to enhance presentation clarity. Specification of nonzero lower bounds is not expected to change the qualitative aspects of subsequent results. The top plot of Figure 3 indicates the energy values used for this example (this historic data corresponds to the PJM Western Hub, Day-Ahead prices for the period of June 1, 2001 through June 10, 2001).22 In the envisioned EMPC implementation, forecasts of energy value will be provided to the controller. In addition to forecast errors, these forecasts are likely to be updated at each time step. However, to avoid the details of forecasting and the closed-loop impacts of forecast error and updates, this example will assume the forecasts to be perfect, in that they are error free and do not change with time. While this assumption is not critical (and will be removed in Section IV), it allows us focus on the fundamental issues and not wonder if observed behavior is due to a forecasting artifact. (As such, in this example, the “hat” notation will be dropped.) The mass of hydrogen in the storage tank is simply the time integral of νs,H2. Then, enforcement of the material balance νs,H2 = νH2 − νG yields the following dynamic model: Ṁ H2 = νH2 −

PG β1

(6)

where β1 = 7.1 MW h/tonne of H2. The discrete-time predictive version of the model is given as M H2(k + 1|i) = M H2(k|i) − bPG(k|i) + c 0 ≤ M H2(k|i) ≤ M Hmax 2 0 ≤ PG(k|i) ≤ PGmax k = i , ..., i + N − 1

(7)

If the sample time is Δts = 1 h, then b = Δts/β1 = 0.141 tonne H2/MW and c = ΔtsνH2 = 90 tonne H2, Mmax H2 = 700 tonne H2 and Pmax G = 1284 MW. It is noted that selecting a value of zero for PG will result in negative values for PN. However, this should not be 3153



Article

Figure 3. Closed-loop simulation of EMPC using a 24-h prediction horizon.

periods in which the energy value is low and the storage is not yet full. Why would the controller perform this way? The answer obviously stems from the abbreviated horizon. As stated in Lima et al.,15 the controller is “myopic”, in the sense that it cannot see beyond its short horizon. As such, the controller decides to use up the hydrogen that it has in inventory in an effort to maximize short-term revenue. However, as illustrated by the 24-h horizon case, if it had known that the value of converting hydrogen to electric energy would soon increase, then it would have continued to stockpile hydrogen rather than use it. One approach to mitigating the inventory creep phenomena is to add a constraint to the final step of the horizon. Specifically, require the predicted inventory at the final step to be greater than some fixed value.15 While this method can reduce the impact of inventory creep, it will also require an adhoc selection of the lower-bound parameters. Also note that, in both the 24-h and 6-h horizon cases, there are intervals in which PG jumps from one extreme to the other and then back to the original extreme. This sort of chattering (or bang−bang-type actuation) seems to occur when the inventory is at an extreme and the energy value is near its average or when the energy value curve experiences sharp changes (for example, at the second, small, daily peak occurring in the early evening). To reduce the occurrence of this chatter, one could add constraints to the first or second time derivative of PG. In Baumrucker and Biegler,14 the square of the first derivative of the manipulated variable is added to the objective function. However, both of these derivative based methods seem to be poorly motivated by the physical situation and will likely compromise optimality of the closed-loop response. Another approach is to force power max output to be in one of three states; 0, Pnom G , or PG . Then, one could add constraints prohibiting two state changes within a certain period (for example, 3 or 4 h). In addition to requiring the use of mixed-integer programming to solve the online EMPC optimization problem, selection of an appropriate waiting period is somewhat arbitrary. Another approach, and the one advocated in this paper, is to extend the horizon size of the EMPC to infinity. However,

an issue, since this power should be available from the grid, and, most likely, at a low price. If Ce(i) denotes energy value during period i, then revenue during period i is Ce(i)PN(i)Δts = Ce(i)PG(i)Δts − Ce(i)(PAC + PCC)Δts. Since the second term cannot be influenced by the manipulated variable, PG(i), an appropriate formulation of the EMPC online optimization is ⎧ Δt ⎨ s M H2(k | i), PG(k | i)⎩ N max

⎪

⎪

i+N−1

∑ k=i

⎫ Ce(k|i)PG(k|i)⎬ ⎭ ⎪

⎪

(8)

subject to M H2(k + 1|i) = M H2(k|i) − bPG(k|i) + c 0 ≤ M H2(k|i) ≤ M Hmax 2 0 ≤ PG(k|i) ≤ PGmax k = i , ..., i + N − 1 M H2(i|i) = M H2(i)

where the constant Δts/N may be dropped, but is retained so that the units of the objective function will be equal to those of average revenue. Clearly, problem (8) can be solved using a standard linear programming (LP) solver. Figure 3 illustrates the closed-loop response of the EMPC policy using a 24-h prediction horizon (N = 24). In this case, the response of the manipulated variable, PG, is as one would desire; at its maximum when energy value is high and at its minimum when energy value is low. The periods in which PG is at the nominal value (Pnom G = 642 MW = β1νH2) correspond to intervals in which the hydrogen storage is full and the energy value is low or the storage is empty and the energy value is high. However, this behavior makes good sense, because the controller has no other reasonable option under these conditions. In the 6-h horizon case (e.g., Figure 4), the EMPC policy seems to lose this good sense. That is, the controller tends to set PG equal to the nominal value more frequently than needed. Specifically, during 3154



Article

Figure 4. Comparison of EPMC using different prediction horizons.

where P̃ G = PG − P̅G, C̃ e = Ce − C̅ e, P̅G = E[PG] and C̅ e = E[Ce]. Equation 10 reflects the general notion that when the energy value is high, power production should be large and production should be small when value is low. The parameter α, yet to be specified, will indicate the magnitude of this relationship. The consequence of eq 10 is that expression (9) can be evaluated as

extending the constrained EMPC problem to an infinite horizon problem would result in an infinite number of variables in the numeric optimization and render the problem intractable. A common approach to similar constrained infinite-time optimal control problems is to enforce only a finite number of constraints.23,24 In this case, the unconstrained infinite-horizon tail is solved analytically and its value function is appended to the constrained finite-time problem as a final cost term, representing the cost-to-go from the final time to infinity. However, if the EMPC objective function for the IGCC problem is applied to an unconstrained problem, the problem will become unbounded and no solution will exist. Thus, our approach will be to modify the unconstrained infinite-horizon tail problem to be statistically constrained in an effort to obtain an analytical solution and approximate the cost-to-go from the final time to infinity.

E[CePG] = E[(Cẽ + Ce̅ )(PG̃ + PG̅ )] = E[Cẽ PG̃ ] + Ce̅ PG̅ = α ΣCe + Ce̅ PG̅

(11)

where ΣCe is the variance of the energy value signal (ΣCe = E[C̃ 2e ]). Thus, if given a characterization of Ce, then one could simply optimize over α and P̅G. (However, for the case being considered, inspection of expression (6) indicates that P̅G must be equal to νH2/β1, regardless of the control policy selected.) A stochastic characterization of Ce begins by recalling the historic energy value data depicted in Figure 3. The data clearly possesses an oscillatory characteristic with a period of 1 day. As such, one could model Ce as the output of an underdamped second-order system driven by white noise (think of a massspring-damper system with a harmonic frequency of 1 day−1). However, experience indicates that such a model will allow too much of the low-frequency energy contained in the white noise to be passed through to Ce. To remove these lowfrequency components, the white noise is first sent through a high pass filter. The net result is the following third-order shaping filter:

III. ECONOMIC LINEAR OPTIMAL CONTROL We now explore the question of developing a linear controller capable of mimicking the EMPC policy. To motivate the discussion, the scenario of Example 1 will be assumed. However, generalization of the approach is straightforward.21 The resulting Economic Linear Optimal Control (ELOC) policy will form the basis of the infinite-horizon EMPC policy to be developed in Section IV. We begin by revisiting the EMPC objective function and consider the impact of taking the limit, with respect to N. ⎧ Δt s lim ⎨ N →∞⎪ N ⎩ ⎪

i+N−1

∑ k=i

⎫ Ce(k|i)PG(k|i)⎬ = E[CePG]Δts ⎪ ⎭ ⎪

ϕ1̇ = ϕ2

(9)

ϕ2̇ = ωc 2(αw1 − ϕ3) − ωc 2ϕ1 − 2χωcϕ2

The result is that the objective function will become the longterm average revenue. The question then becomes this: Can one calculate this long-term average, and can one select a policy for PG such that this average is maximized? One approach is to assume PG has the following form:

PG̃ = αaCẽ

ϕ3̇ = Ce =

(10) 3155

αw1 − ϕ3 τh ϕ1 α

+ Ce̅

(12)



Article

Figure 5. Comparison of the ELOC policy and the EMPC policy (24-h prediction horizon).

where ωc = 2π/τc, τc = 24 h, χ = 0.1, τh = 1 h, and w1 is a Gaussian, zero-mean white noise process with the following spectral density: ⎛ 4χ ⎞⎛ ω 2τ 2 + 2χωcτh + 1 ⎞ ⎟ΣCe Sw1 = ⎜ ⎟⎜ c h ωc 2τh 2 ⎝ ωc ⎠⎝ ⎠

where M̅ H2 is arbitrarily defined as MHmax /2. The process 2 (augmented with the shaping filter) then can be written as x ̇ = Ax + Bu + αGw z = Dx x + Duu

(13)

z min ≤ z ≤ z max

The idea behind this shaping filter model is that, for all parameter values (τc, χ, τh, C̅ e, and α), a stochastic covariance analysis on expression (12) will indicate that the calculated variance of Ce is equal to the ΣCe parameter value selected in expression (13). As such, it is highlighted that the τc, χ, τh values selected above are not special, but simply represent parameter values such that realizations of Ce match, purely on a visual basis, with the data of Figure 3. While this is sufficient for our purpose of illustrating the method, a more rigorous system identification approach is likely warranted, which may include changing the model structure from that given in expression (12). The top plot of Figure 5 is a realization of Ce generated from expression (12), using the above parameters, as well as α = 1 MW2 h/$, C̅ e = 20 $/(MW h), and ΣCe = 12 $/(MW h)2. An important aspect of the definition of Ce is the fact that ϕ1 is equal to αC̃ e. Thus, to enforce an approximation of condition (10), one could require

E[(PG̃ − ϕ1)2 ] < ε

(16)

where ⎡0 0 0 0 ⎤ ⎥ ⎢ 0 1 0 ⎥ ⎢0 A=⎢ 0 − ωc 2 − 2χωc − ωc 2 ⎥ ⎥ ⎢ ⎢⎣ 0 0 0 − 1/τh ⎥⎦ ⎡ 0 ⎤ ⎥ ⎢ ⎢ 0 ⎥ G=⎢ 2⎥ ω ⎢ c ⎥ ⎢⎣1/τh ⎥⎦

z max

where ε is selected to be sufficiently small, but not so small that numeric feasibility becomes an issue. However, since PG is the manipulated variable, satisfaction of this inequality should be achievable for rather small values of ε (ε is set to 0.01 MW2 for all cases). The next step is to recast the compound system into deviation variable form:

(17)

⎡0⎤ ⎢ ⎥ Du = ⎢ 1 ⎥ ⎣1 ⎦

z min

⎡ M Hmax − M̅ H ⎤ 2 2 ⎢ ⎥ = ⎢ P max − P ̅ ⎥ G ⎢ G ⎥ ⎣ ⎦ ε

⎡− M̅ H ⎤ 2 ⎢ ⎥ = ⎢ − P̅ ⎥ G ⎢ ⎥ ⎣ −ε ⎦

(18)

The continuous-time model can then be converted to a discretetime form, using the sample and hold method (Ad = eAΔts, 25 Δts At s At Bd = ∫ Δt 0 e B dt, Gd = ∫ 0 e G dt, and Σw = Sw/Δts):

T x = ⎡⎣ M H2 − M̅ H2 ϕ1 ϕ2 ϕ3 ⎤⎦

x(i + 1) = Ad x(i) + Bd u(i) + αGdw(i)

u = [PG − PG̅ ]

z(i) = Dx x(i) + Duu(i)

w = [w1]

z = ⎡⎣ M H2 − M̅ H2 PG − PG̅ (PG − PG̅ ) − ϕ1⎤⎦T

1/β1⎤ ⎥ 0 ⎥ ⎥ 0 ⎥ 0 ⎦

Sw = [Sw1]

⎡1 0 0 0⎤ ⎢ ⎥ Dx = ⎢ 0 0 0 0⎥ ⎣0 − 1 0 0⎦

(14)

⎡− ⎢ ⎢ B=⎢ ⎢ ⎣

z min ≤ z(i) ≤ z max

(15) 3156

(19)



Article

The ELOC controller is then assumed to be a simple linear feedback of the state: u(i) = Lx(i). Given a candidate controller (L), the variance of the jth performance variable (ζj) is calculated as ζj =

ρj ΣzρjT

Application of Theorem 1 to eq 24 results in the following convex optimization: max

Y , X ≥ 0, μj , α

(20) T

Σz = (Dx + DuL)Σx(Dx + Dx L)

(21)

Σx = (Ad + Bd L)Σx(Ad + Bd L)T + α 2Gd ΣwGd T

(22)

ζj < ( −z jmin)2

and

max

⎡ μj ρj (Dx X + DuY )⎤ ⎢ ⎥ ⎢ ⎥>0 T T + ρ ( ) D X D Y X u j ⎣ x ⎦ μj < zj̅ 2

(23)

{αΣC + Ce̅ PG̅ } e

zj̅ = min{z jmax , −z jmin},

(24)

Σx = (Ad + Bd L)Σx(Ad + Bd L)T + α 2Gd ΣwGd T ζj = ρj (Dx + DuL)Σx(Dx + Dx L)T ρj T < zj̅ 2 j = 1 , ..., 3

To arrive at a convex formulation of expression (24), the following theorem (a slight generalization of Theorem 6.1 from Chmielewski and Manthanwar26) is applied. Theorem 1: There exists stabilizing L, Σx ≥ 0, and ζj ≥ 0, j = 1 , ..., nz, such that Σx=(Ad + Bd L)Σx(Ad + Bd L)T + α 2Gd ΣwGd T

ζj = ρj (Dx + DuL) ∑ (Dx + DuL)T ρj T ,

(25)

L* = [ 0.067 − 0.018 0.017 − 0.0009 ] j = 1 , ..., nz (26)

ζj
0, T T ( + ) ρ D X D Y X x u j ⎣ ⎦

(28)

j = 1 , ..., nz (29)

and μj < zj̅ 2 ,

j = 1 , ..., nz

(33)

The second two plots of Figure 5 compare the resulting ELOC policy with the EMPC policy using the energy value data of the first plot of Figure 5. As expected the ELOC policy does not enforce constraints pointwise-in-time, but does enforce them statistically. However, the important point is that the ELOC trajectory has a general shape that is remarkably similar to the EMPC trajectory. Note that a trajectory with less pointwise-in-time constraint violations could be achieved by tightening the statistical constraints of the ELOC. (For example, in expression (31), μj < zj̅ 2 could be replaced with μj < (zj̅ /2)2.) This, however, would still not guarantee that pointwise-in-time constraints would not be violated; it only ensures that they would be violated less frequently. Furthermore, if the statistical constraints are tightened too much, then the trajectory would be significantly altered from that of the EMPC. The next section will illustrate a novel approach to enforce pointwise-in-time constraints while preserving the positive features of the ELOC.

if and only if there exists X > 0, Y and μj > 0, j = 1 , ..., nz, such that ⎡ X (Ad X + Bd Y ) αGd ⎤ ⎢ ⎥ ⎢ ⎥ ⎢(A X + B Y )T X 0 ⎥>0 d ⎢ d ⎥ ⎢ ⎥ ⎢ αGd T Σw−1⎥⎦ 0 ⎣

(32)

After rescaling, with respect to α*, the ELOC feedback gain is determined to be

x

z ̅j2 ,

j = 1 , ..., 3

Because of the linearity of the objective function and the convexity of the constraints, any local optimum of problem (31) is guaranteed to be a global optimum. It is additionally noted that if the solution to expression (31) is X* and Y*, then L* = Y*(X*)−1 is the solution to expression (24), within the accuracy of the strict inequality constraints of expression (31). To arrive at the ELOC policy, the solution L* must be scaled, with respect to α*. That is, the columns of L* corresponding to the shaping filter states (the second through fourth states) must be multiplied by α*. This rescaling will make the policy appropriate for a price of electricity generated by the equations described by expression (12), with α set equal to 1. Example 2a: Reconsider the scenario of Example 1, and assume that the shaping filter used to model Ce has parameter values as indicated in this section. Problem (31) is solved using the Linear Matrix Inequality (LMI) optimization routines found in the Robust Control Toolbox of Matlab. The solution to problem (31) is found to be α* = 774 MW2 h/$, and

subject to

zj̅ = min{zjmax , − zjmin},

(31)

⎡ X (Ad X + Bd Y ) αGd ⎤ ⎢ ⎥ ⎢(A X + B Y )T ⎥>0 X 0 d ⎢ d ⎥ ⎢⎣ αGd T Σw−1⎥⎦ 0

The next step is to find a linear feedback, L, to maximize our expression for average revenue (see eq 11). The following nonlinear optimization is posed to determine such a feedback. L , Σx ≥ 0, ζj , α

e

subject to

where ρj is the jth row of an identity matrix and Σx is the positive definite solution of eq 22. Rather than enforce pointwise-in-time constraints, zmin ≤ zj(i) ≤ zmax j j , the central aspect of the ELOC method is to follow a pseudo- (or chance-) constrained approach by enforcing the following statistical constraints: ζj < (z jmax )2

{αΣC + Ce̅ PG̅ }

(30) 3157



Article

Figure 6. Comparison of the IH-EMPC and the EMPC policies using a 24-h horizon.

IV. INFINITE-HORIZON ECONOMIC MPC We now turn to the question of converting the ELOC policy into one that is capable of enforcing constraints pointwise-in-time. To do so, we revisit the quadratic version of MPC and highlight the following important property. Consider an unconstrained finitehorizon MPC controller:

The second important point is that, for some linear feedback, there exists an equivalent unconstrained MPC controller. Thus, if such a property holds for the feedback generated by expression (31), then the addition of constraints to its MPC counterpart will generate a policy similar to the original linear feedback. In the earlier work by Chmielewski and Manthanwar,26 it is shown that a class of problems (to which expression (31) is a member) is guaranteed to possess an equivalent unconstrained MPC controller. Furthermore, Chmielewski and Manthanwar26 provided the following theorem on inverse optimality, which can be used to synthesize weighting matrices Q, R, M, and P (from L, Ad, and Bd), such that the resulting unconstrained MPC is equivalent to the original policy u(i) = Lx(i). Theorem 6.2 of Chmielewski and Manthanwar26 is restated here, for the sake of convenience. Theorem 2: If there exists P > 0 and R > 0, such that

T ⎧i+N−1 ⎡ x(k|i)⎤ ⎡ Q M ⎤ ⎪ ⎥ ⎢ ⎥[x(k|i)u(k|i)] min ⎨ ∑ ⎢ T x(k | i), u(k | i)⎪ ⎢ ⎥ ⎩ k = i ⎣ u(k|i)⎦ ⎣ M R ⎦

⎫ ⎪ + x(i + N |i) Px(i + N |i))⎬ ⎪ ⎭ T

(34)

subject to x(k + 1|i) = Ad x(k|i) + Bd u(k|i),

k = i , ..., i + N − 1

⎡ P − Ad TPAd + LT LT(R + Bd T PBd ) + Ad TPBd ⎤⎥ ⎢ T (R + Bd PBd )L ⎢ ⎥>0 ⎢ ⎥ T T ⎢⎣(R + Bd PBd )L + Bd PAd ⎥⎦ R

x(i|i) ≤ x(i)

The property then is stated as follows: If P is selected such that P = Ad TPAd + Q

(37)

− (M + Ad TPBd )(R + Bd T PBd )−1(M + Ad TPBd )T

then Q =̂ P − ATd PAd + LT(R + BTd PBd)L and M =̂ LT(R + BTd PBd) + ATd PBd will be such that

(35)

then the MPC controller of expression (34) will be independent of horizon size (N) and will be identical to using the following linear feedback gain: L = −(R + Bd T PBd )−1(M + Ad TPBd )T

⎡ Q M⎤ ⎢ ⎥>0 ⎣ MT R ⎦

(36)

and P and L satisfy eqs 35 and 36. In summary, construction of the infinite-horizon EMPC is as follows. Begin by solving problem (31) to determine LELOC. Then, implement the procedure suggested by Theorem 2 to synthesize weighting matrices Q, R, M, and P, which are companions to LELOC. Finally, implement the following constrained MPC policy using a standard QP solver.

An equivalent statement would be: u*(i|i) = Lx(i). The first important point is the independence of the policy, with respect to horizon size, which is made possible by the appropriate selection of P. In essence, selecting P based on eq 35 makes all unconstrained finite-horizon MPC policies equal to the infinite-horizon policy, which is the classic linear quadratic regulator feedback given in eq 36. 3158



Article

Figure 7. Comparison of the IH-EMPC policies using different horizons. T ⎧i+N−1 ⎡ x(k|i)⎤ ⎡ Q M ⎤ ⎪ ⎥ ⎢ ⎥[x(k|i)u(k|i)] min ⎨ ∑ ⎢ T x(k | i), u(k | i)⎪ ⎥ ⎢ ⎩ k = i ⎣ u(k|i)⎦ ⎣ M R ⎦ ⎫ ⎪ + x(i + N |i)T Px(i + N |i))⎬ ⎪ ⎭

different weights. However, the important aspect is to verify that a substitution of the weights back into eqs 35 and 36 regenerates the policy LELOC. Using a 24-h horizon and the energy values depicted in the top plot of Figure 6 (the same as those given in Figure 5), the EMPC policy is compared to the infinite-horizon EMPC (IH-EMPC). Clearly, there is a strong resemblance between the two policies. More remarkable is the fact that the chatter observed in the EMPC seems be damped out almost completely. That is, the IH-EMPC policy shows a tendency toward the chatter actions but is slowed by the move suppression aspects of the quadratic objective. The plots of Figure 7 illustrate the policy’s insensitivity, with respect to horizon size. The fact that even with a horizon of 1 h (one time-step of prediction), the trajectory is virtually unchanged indicates the strong influence of the final cost term: x(i + N|i)TPx(i + N|i)).

(38)

subject to x(k + 1|i) = Ad x(k|i) + Bd u(k|i) + Gdŵ(k|i), k = i , ..., i + N − 1 z(k|i) = Dx x(k|i) + Duu(k|i), z min ≤ z(k|i) ≤ z max ,

k = i , ..., i + N − 1

k = i , ..., i + N − 1

x(i|i) ≤ x(̂ i)

Example 2b: Continuing Example 2a, the linear feedback described by expression (33) generates the weighting matrices: ⎡ 15.7 − 3304 3021 − 162 ⎤ ⎢ ⎥ ⎢ ⎥ 5 5 − 3304 6.96 × 10 − 6.37 × 10 34200 ⎢ ⎥ Q=⎢ ⎥ ⎢ 3021 − 6.37 × 105 5.83 × 105 − 31300 ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ − 162 34200 − 31300 1684 ⎦

V. PARTIAL STATE INFORMATION IH-EMPC While the plots of Figures 6 and 7 are impressive, the scenario of Example 2 provides the IH-EMPC policy with two unrealistic advantages; perfect forecasting and the generation of energy value data by the shaping filter model contained within the IH-EMPC policy. The assumption of a perfect forecast boils down to assuming perfect knowledge of the white noise sequence w(k) for k > i, where i is the current time. In terms of problem (38), that is ŵ (k|i) = w(k) for all i. In the imperfect forecasting case, the procedure is as follows. Using a shaping filter model, such as that in expression (12) with α = 1, apply a state estimator to determine x̂(i). Then, implement a state predictor (essentially a simulation of expression (12) with x(i|i) = x̂(i) and ŵ (k|i) = 0) to generate forecasts x(k|i). Application of this procedure to the historic data of Example 1 (top plot of Figure 3) resulted in the plots of Figure 8. Clearly, both policies fail. The problem seems to be an overcompensation for the short-term average energy value. That is, the controllers hoard H2 in the first half of the simulation (t < 10 days), because of the fact that Ce is less than C̅ e (= 20 $/(MW h)), and deplete the H2 inventory in the second half, when Ce is mostly greater than C̅ e. One way to address this issue is to restructure the shaping filter to capture the changes in the short-term average. For example, the filter described in expression 12 could be augmented with a

R = [3482]

(39)

⎡ − 234 ⎤ ⎢ ⎥ ⎢ 49300 ⎥ ⎢ ⎥ M=⎢ ⎥ ⎢− 45000 ⎥ ⎢ ⎥ ⎢⎣ 2418 ⎥⎦ ⎡ 7.26 × 10−2 2.2 12.5 0.2 ⎤ ⎢ ⎥ 2.2 498 540 − 19 ⎥ ⎢ × 10−3 P=⎢ 12.5 540 8760 − 224 ⎥ ⎢ ⎥ ⎣ − 19 − 224 6660 ⎦ 0.2 (40)

Note that the weighting matrices Q, R, M, and P are not unique, in that identical implementations of Theorem 2 on different computers or different versions of Matlab will likely result in 3159



Article

Figure 8. Application of historic energy value data to the EMPC and IH-EMPC policies using expression (12) as the forecasting model.

low-pass filter driven by an independent white noise process. In this case, expression 12 would be replaced with

short-term average. However, in reality, the state estimator cannot perform this separation perfectly. Specifically, the estimator has available to it only the energy value curve, which it models as the sum of ϕ1 and ϕ4. Then, it must separate the two based only on frequency content. Thus, the problem with our previous IH-EMPC policy is that the quality of the state estimator has not been communicated to the design procedure. To be able to communicate this estimator quality information, a Partial State Information (PSI) formulation of the ELOC policy must be developed. We begin with a review of state estimation. If given a process model x(i + 1) = Adx(i) + Bdu(i) + Gdw(i) and a measurement equation y(i) = Cx(i) + Fv(i), where w and v are zero mean, Gaussian white noise sequences, then the optimal (minimum error variance) state estimate, x̂, is generated from the following system:28

ϕ1̇ = ϕ2 ϕ2̇ = ωc 2(α1w1 − ϕ3) − ωc 2ϕ1 − 2χωcϕ2 ϕ3̇ = ϕ4̇ = Ce =

(α1w1 − ϕ3) τh α2w2 − ϕ4 τl ϕ1 α1

+

ϕ4 α2

+ Ce̅

(41)

where τl = 6 h and w2 is a Gaussian, zero-mean white noise process that is independent of w1 and the spectral densities of the two are defined as follows, with γ = 0.85: ⎛ 4χ ⎞⎛ ω 2τ 2 + 2χωcτh + 1 ⎞ ⎟ΣCe Sw1 = γ ⎜ ⎟⎜ c h ωc 2τh 2 ⎝ ωc ⎠⎝ ⎠

(42)

Sw2 = (1 − γ )2τl ΣCe

(43)

x(̂ i + 1) = Ad x(̂ i) + Bd u(i) + K (y(i) − Cx(̂ i))

(44)

where K = Ad ΣeCT(C ΣeCT + F ΣvFT)−1

One could then define z3 as P̃G − ϕ1 − ϕ4 and let the optimization select α1 and α2 such that the resulting policy will give appropriate weight to each state of the shaping filter (see Mendoza-Serrano and Chmielewski27 for details on the implementation of this procedure). The result of the optimization is as one would expect: α2* is found to be zero, indicating that the ELOC will place all its emphasis on short-term variability and ignore the short-term average. The result of using eq 41 as the forecasting model is shown in Figure 9. In this case, both policies perform much better. However, both still seem to have trouble when the short-term average is significantly above or below C̅ e. That is, in the first three days, the policies hold on to the H2 inventory for too long, and during the periods of 10−15 days and 17−20 days, the policies do not store enough H2. The problem stems from the fact that the policies think that the short-term variability can be perfectly distinguished from the

(45)

and Σe is the positive definite solution to Σe = Ad ΣeAdT − Ad ΣeCT(C ΣeCT + F ΣvFT)−1C ΣeAdT + Gd ΣwGdT

(46)

The following, easily verified, facts will be utilized in the subsequent development. Fact 4.1: If the matrices Ad, C, Gd, Σw, F, and Σv are block diagonal (with appropriate block dimensions), then the matrices K and Σe will be block diagonal. Specifically, if Ad = diag{Adj}, C = diag{Cj}, Gd = diag{Gdj}, Σw = diag{Σwj}, F = diag{Fj}, and Σv = diag{Σvj}, then K = diag{Kj} and Σe = diag{Σej}, where Kj = Adj ΣejCj T(Cj ΣejCj T + Fj ΣvjFj T)−1 3160

(47)



Article

Figure 9. Application of historic energy value data to the EMPC and IH-EMPC policies using eq 41 as the forecasting model.

Thus, given the estimator system

Σej = Adj ΣejAdjT − Adj ΣejCj T(Cj ΣejCj T + Fj ΣvjFj T)−1Cj ΣejAdj T + Gdj ΣwjGdj T

x(̂ i + 1) = (Ad + Bd L)x(̂ i) + η(i)

(48)

z(i) = Dx x(i) + DuLx(̂ i) + Dx ̂x(̂ i)

Fact 4.2: If the matrices Ad, C, Gd, and F are defined as ⎡ Ad 0 0 ⎤ ⎥, Ad = ⎢ ⎢⎣ 0 A ⎥⎦ d0

the variance of the jth element of z is calculated as

⎡C 0 0 ⎤ ⎥, C=⎢ ⎢⎣ 0 C0 ⎥⎦

⎡ α1Gd 0 ⎤ ⎥ , and Gd = ⎢ ⎢⎣ α2Gd 0 ⎥⎦

⎡ α1F0 ⎤ ⎥ F=⎢ ⎣ α2F0 ⎦

ζj = ρj ΣzρjT

(54)

Σz = (Dx + DuL + Dx )̂ Σ x(̂ Dx + DuL + Dx )̂ T + Dx ΣeDxT (55)

(49)

T

Σ x ̂ = (Ad + Bd L)Σ x(̂ Ad + Bd L) + Ση

then ⎡ α 2Σ ⎤ 1 e 0 α1α2 Σe 0 ⎥ Σe = ⎢ 2 ⎢α α Σ ⎥ ⎣ 1 2 e0 α2 Σe0 ⎦

and

K = [ α1K 0 α2K 0 ]

Σ x ̂ = (Ad + Bd L)Σ x(̂ Ad + Bd L)T + Ad ΣeCT(C ΣeCT + Σv)−1C ΣeAd T

where K0 and Σe0 are from eqs 47 and 48, with j = 0. Furthermore,

× ⎣⎡α1C0 Σe0Ad 0

α2C0 Σe0Ad 0

(58)

and ζj < z ̅j2 , j = 1 , ..., nz

(59)

if and only if there exists X > 0, Y and μj > 0, j = 1 , ..., nz, such that

(51)

⎡ (Ad X + Bd Y ) (Ad ΣeCT) ⎤ X ⎢ ⎥ ⎢(A X + B Y )T ⎥>0 0 X d ⎢ d ⎥ ⎢ (A Σ CT)T 0 (C ΣeCT + Σv)⎥⎦ ⎣ d e

Covariance analysis under the PSI feedback structure, u(i) = Lx̂(i), utilizes two important properties of optimal state estimation. First, the optimal estimate, x̂(i), and the error signal, e(i) =̂ x(i) − x̂(i), are orthogonal, in that E[x̂(i)e(i)] = 0. An immediate consequence of this fact is that Σx̂ = Σx − Σe. The second important property is that the innovations sequence, defined as η(i) =̂ K(y(i) − Cx̂(i)), is a zero mean, white noise process with covariance Ση = Ad ΣeCT(C ΣeCT + F ΣvF )−1C ΣeAd T

(57)

ζj = ρj (Dx + DuL + Dx )̂ Σ x(̂ Dx + DuL + D x )̂ T ρj T + ρj Dx ΣeDxTρj T

Ad ΣeCT(C ΣeCT + F ΣvFT)−1C ΣeAd T ⎡α A Σ C T⎤ 1 d 0 e0 0 ⎥ (C Σ C T + F0 ΣvF0 T)−1 =⎢ ⎢ α A Σ C T ⎥ 0 e0 0 ⎣ 2 d 0 e0 0 ⎦ T⎤ ⎦

(56)

Theorem 3: There exists stabilizing L, Σx̂ ≥ 0 and ζj ≥ 0, j = 1 , ..., nz, such that (50)

T

(53)

(60)

⎡ μj − ρj Dx ΣeDxTρj T ρj ((Dx + Dx )̂ X + DuY )⎤ ⎢ ⎥ ⎢ ⎥ > 0, T T X ⎢⎣((Dx + Dx )̂ X + DuY ) ρj ⎥⎦ j = 1 , ..., nz

(52) 3161

(61)



Article

The second point concerns the use of α1 and α2 in expressions (63) and (64), which will allow the optimization to select how much weight it will give to the short-term variability, ϕ̂ 1/α1, and how much to give to the short-term average, ϕ̂ 4/α2. An evaluation of the average revenue in the context of expression (67), yields the following expression:

and μj < z ̅j2 ,

j = 1 , ..., nz

(62)

Proof: If the matrices Σx, αGd, Σw, Dx, and ζi of Theorem 1 are redefined as ∑x̂, AdΣcCT, (CΣcCT + Σv)−1, Dx + Dx̂, and ζi − ρiDxΣxDTx ρTi , then Theorem 3 will result. ■ Example 3: This example will continue the scenarios of Examples 1 and 2, but update the shaping filter model to eqs 41 and extend the ELOC method to the PSI framework. Consider the following block diagonal process model: ⎡ Ad 0 0 ⎤ ⎢ 0 ⎥ Ad = ⎢ 0 A d1 0 ⎥ ⎢ ⎥ 0 A d1 ⎥⎦ ⎢⎣ 0

⎡ Bd ⎤ ⎢ 0⎥ Bd = ⎢ 0 ⎥ ⎢⎣ 0 ⎥⎦

E[PG̃ Cẽ ] = d1α1 + d 2α2

where d1 = ρ1Σx̂1ρT1 + ρ1Σx̂1ρT4 , d2 = ρ1Σx̂1ρT4 + ρ4Σx̂1ρT4 , Σx̂1 = Σx1 − Σe1 (Σe1 is the positive definite solution to eq 48, with j = 1, and Σx1 satisfies Σx1 = Ad1Σx1ATd1 + Gd1Σw1GTd1. Finally, using Theorem 3, along with Facts 4.1 and 4.2, the PSI version of the ELOC design problem, in convex form, is stated as

⎡G d 0 ⎤ ⎢ 0 ⎥ Gd = ⎢ 0 α1Gd1 ⎥ ⎢ ⎥ ⎢⎣ 0 α2Gd1 ⎥⎦

⎡ Σw 0 ⎤ 0 ⎥ Σw = ⎢ ⎢⎣ 0 Σw ⎥⎦ 1

max

{d1α1 + d 2α2}

Y , X ≥ 0, μj , α1, α2

⎡ F0 0 ⎤ ⎢ ⎥ F = ⎢ 0 α1F1 ⎥ ⎢ ⎥ ⎣ 0 α2F1⎦

⎡ (Ad X + Bd Y ) M1 ⎤ X ⎢ ⎥ ⎢(A X + B Y )T 0 ⎥>0 X d d ⎢ ⎥ ⎢⎣ 0 M1T M 0 ⎥⎦

⎡ Σv 0 ⎤ 0 ⎥ Σv = ⎢ ⎢⎣ 0 Σv ⎥⎦ 1

⎡ μ − ρ D Σ DTρ T ρj ((Dx + Dx )̂ X + DuY )⎤ j j x e x j ⎢ ⎥ ⎢ ⎥>0 T X ⎢⎣((Dx + Dx )̂ X + DuY ) ρj ⎥⎦

(64)

where the elements of the first block are Ad0 = 1, Bd0 = −Δts/β1, Gd0 = Δts, C0 = F0 = 1, Σw0 = Σv0 = (0.1)2/Δts; and the elements of s A1t the second block are Ad1 = eA1Δts, Gd1 = ∫ Δt 0 e G1 dt, Σw1 = Sw/Δts, and Σv1 = Sv1/Δts, where ⎡ 0 1 0 ⎢ 2 2 ⎢− ωc − 2χωc − ωc A1 = ⎢ − 1/τh 0 ⎢ 0 ⎢ − 0 0 ⎣ 0

(69)

subject to

(63)

⎡C 0 0 0 ⎤ ⎢ ⎥ ⎢ C = 0 C1 0 ⎥ ⎢ ⎥ ⎢0 0 C⎥ ⎣ 1⎦

⎤ ⎡(C Σ C T + Σ ) 0 v0 ⎥ ⎢ 0 e0 0 M0 = ⎢ ⎥ T 0 ( ) C C Σ + Σ v1 ⎦ 1 e1 1 ⎣

⎤ ⎥ 0 ⎥ ⎥ 0 ⎥ ⎥ 1/τl ⎦ 0

⎤ ⎡A Σ C T 0 ⎥ ⎢ d0 e0 0 ⎥ ⎢ 0 α1A d1Σe1C1T ⎥ M1 = ⎢ ⎥ ⎢ 0 α2A d1Σe1C1T ⎥⎦ ⎢⎣ μj < zj̅ 2

⎡ 0 0 ⎤ ⎢ 2 ⎥ 0 ⎥ ⎢ ωc G1 = ⎢ ⎥ ⎢1/τh 0 ⎥ ⎢ ⎥ ⎣ 0 1/τl ⎦

zj̅ = min{z jmax , − z jmin},

⎡ Sw1 0 ⎤ ⎥ Sw = ⎢ ⎢⎣ 0 Sw2 ⎥⎦

(65)

⎡ 0 0 0 0⎤ ⎢ ⎥ Dx1̂ = ⎢ 0 0 0 0 ⎥ ⎣− 1 0 0 0 ⎦

and

⎡0 0 0 0⎤ ⎢ ⎥ Dx2̂ = ⎢ 0 0 0 0⎥ ⎣ 0 0 0 − 1⎦

L ELOC = [ 0.0445 86.58 − 7.861 0.5385 0 0 0 0 − 45.13]

(70)

The fact that α2* turns out to be negative indicates that the policy recognizes a correlation between ϕ̂ 1 and ϕ̂ 4, in the sense that when the short-term average of ϕ̂ 1 is greater than zero, the short-term average of ϕ̂ 4 will also be positive. Thus, the selection of a negative α2*, which manifests as a negative value for the ninth element of LELOC, has the effect of canceling out the short-term average of ϕ̂ 1. The plots of Figure 10 compare the PSI and FSI versions of the IH-EMPC policy. While these trajectories seem to be very similar, the PSI policy uses more of the dispatch capabilities. To see this improvement, the closed-loop revenue using 3 months of historic data (June 2001−August 2001, from PJM22) was calculated. As seen in Table 1, the revenue increase resulting from the PSI policy is 2.1 percentage points greater than that of the FSI policy (24.0%, vs 21.9%). It is also worth noting that the PSI policy achieves almost 90% of the revenue increase that would be realized if given perfect forecasting.

(66)

min

Finally, the z and z vectors are as described in Example 2. The third row of the performance output indicates an enforcement of the following constraint: E⎡⎣(PG̃ − ϕ1̂ − ϕ4̂ )2 ⎤⎦ < ε 2

j = 1 , ..., 3

The solution to problem (69) is found to be α1* = 87.3 MW2 h/$ and α2* = −47.0 MW2 h/$ and

Sw1 and Sw2 are as described in eqs 42 and 43, C1 = [1 0 0 1], F1 = 1, and Sv1 = (0.1)2. The performance output is defined as Dx = [Dx0 Dx1 Dx2], Du = [0 1 1]T, Dx̂ = [Dx̂0 Dx̂1 Dx̂2], where Dx0 = [1 0 0]T, Dx1 = Dx2 = Dx̂0 = 0,

max

(68)

(67)

The first point about expression (67) is that the manipulated variable, PG, is being forced to track the estimates of the shaping filter states. This is appropriate, since the controller does not have access to the “true” state values, which, in reality, do not even exist. 3162



Article

Figure 10. Comparison of PSI and FSI versions of the IH-EMPC policy using eq 41 as the forecasting model.

Table 1. Calculated Policy Revenues over a Three-Month Period (All Using a 24-h Horizon) case

revenue (× 103 $/day)

revenue increase (%)

no dispatch EMPC: perfect forecast EMPC: imperfect forecast IH-EMPC: imperfect forecast (FSI) IH-EMPC: imperfect forecast (PSI)

557.6 765.7 683.9 713.5 736.0

27.2 18.5 21.9 24.0

disturbances other than the value of electricity. In MendozaSeriano and Chmielewski,12,27 a relation similar to expression (10) included a second weighted term that was proportional to a signal orthogonal to C̃ e. This additional degree of freedom allowed the manipulated variable to range from full correlation with C̃ e to zero correlation with C̃ e. Another important question concerns the assumption of the ELOC being restricted to an unstructured linear feedback of the state, or state estimate. Generally, this assumption is required to guarantee that problem (31) will be convex. For example, if a nonlinear controller was allowed, then eqs 20−22 would have to be replaced with Fokker−Planck relations that could not be converted to convex constraints. It is also noted that restricting the linear feedback gain (L) to being, in some way, structured will again preclude the conversion of eqs 20−22 to a convex form. However, as highlighted in Section V, suitable selection of the shaping filter can provide an approximation of integral action. In Section V, we jumped to the case of imperfect forecasting to illustrate how to implement the IH-EMPC using energy value data that are not generated by the shaping filter. Thus, the question remains regarding how to implement the IH-EMPC if perfect forecasts are available (or, alternatively, if high-quality forecasts from a black-box model, unavailable to the IH-EMPC, are provided). One approach is to extend the methods of Section V. Specifically, convert the state predictor to a state smoother by considering forecasts as measurements (at future times) to obtain a higher-quality sequence x(k|i). In this case, the variance of the noise added to future measurements could be set to commensurate with their quality (small for the near future and large for the distant future). The downside of this approach is that the shaping filter may not have enough structure to capture the details of the high-quality forecasts, especially if the variance of the measurement noise is assumed very small, and will result in suboptimality due to modifications of the forecasts. Thus, an alternative approach is to modify problem (38) to a form similar to that of expression (9), but with the terminal cost term included. In this case, the electricity value used in the terminal cost term

■

CONCLUSIONS In this work, we have introduced the notion of EMPC as a method to maximize revenue from a dispatch-capable IGCC power plant. However, this EMPC policy was found to possess several shortcomings, including inventory creep and the existence of bang−bang-type actuation. To alleviate these issues, a new economic-based linear optimal control policy (ELOC) was developed, and through the use of inverse optimality on the ELOC policy, an infinite-horizon formulation of the EMPC policy was developed. This policy was shown to possess very attractive characteristics, but seemed to fail if the energy value data possessed nonzero short-term averages. To address this issue, a partial state information version of the IH-EMPC policy was developed and shown to perform quite well in the face of nonzero short-term averages. Concerning the ELOC policy, it is highlighted that expression (10) is critical to achieving a convex approach to the determination of the ELOC feedback gain. However, as illustrated in Section V, a modification of this relation, from expression (10) to expression (67), was required to address nonzero short-term average energy values, as well as the PSI framework. What we can say with certainty is that the imposition of relations similar to expression (10) will be case-specific and it is an open question as to the existence of a generic relation applicable to all situations. Consider, for example, the case of a system being driven by 3163



Article

(14) Baumrucker, B. T.; Biegler, L. T. MPEC strategies for cost optimization of pipeline operations. Comput. Chem. Eng. 2010, 34, 900− 913. (15) Lima, R. M.; Grossmann, I. E.; Jiao, Y. Long-term scheduling of a single-unit multi-product continuous process to manufacture high performance glass. Comput. Chem. Eng. 2011, 35 (3), 554−574. (16) Kostina, A. M.; Guillén-Gosálbeza, G.; Meleb, F. D.; Bagajewiczc, M. J.; Jiméneza, L. A novel rolling horizon strategy for the strategic planning of supply chains: Application to the sugar cane industry of Argentina. Comput. Chem. Eng. 2011, 35, 2540−2563. (17) Diehl, M.; Amrit, R.; Rawlings, J. B. A Lyapunov function economic optimizing model predictive control. IEEE Trans. Autom. Control 2011, 56 (3), 703−707. (18) Huang, R.; Biegler, L. Stability of economically-oriented NMPC with periodic constraint. In Preprints of the 18th IFAC World Congress, Milano, Italy, 2011; pp 10499−10504. (19) Heidarinejad, M.; Liu, J.; Christofides, P. D. Economic model predictive control of nonlinear process systems using Lyapunov techniques. AIChE J. 2012, 58, 855−870. (20) Klara, J. Cost and performance baseline for fossil energy plants, Volume 1: Bituminous coal and natural gas to electricity. Technical Report DOE/NETL-2007/1281; National Energy Technology Laboratory (NETL): Washington, DC, 2007. (21) Yang, M. W.; Omell, B. P.; Chmielewski, D. J. Controller design for dispatch of IGCC power plants. In Proceedings of the American Control Conference, Montreal, Canada, June 27−29, 2012; pp 4281− 4286. (22) PJM. PJM Day-Ahead Historical Data. http://www.pjm.com/ markets-and-operations/energy/day-ahead/day-ahead-historical.aspx (accessed on March 30, 2012). (23) Chmielewski, D. J.; Manousiouthakis, V. On constrained infinitetime linear quadratic optimal control. Syst. Control Lett. 1996, 29, 121. (24) Manousiouthakis, V.; Chmielewski, D. J. On constrained infinitetime nonlinear optimal control. Chem. Eng. Sci. 2002, 57, 105−114. (25) Burl, J. B. Linear Optimal Control; Addison−Wesley: Menlo Park, CA, 1999. (26) Chmielewski, D. J.; Manthanwar, A. M. On the tuning of predictive controllers: Inverse optimality and the minimum variance covariance constrained control problem. Ind. Eng. Chem. Res. 2004, 43, 7807−7814. (27) Mendoza-Serrano, D. I.; Chmielewski, D. J. Controller and system design for HVAC with Thermal Energy Storage. In Proceedings of the American Control Conference, Montreal, Canada, June 27−29, 2012; pp 3669−3674. (28) Anderson, B. D. O.; Moore, J. B. Optimal Filtering; Prentice Hall: Englewood Cliffs, NJ, 1979.

would need to be the state estimate at time i + N, based on the forecast provided up to time i + N. This type of controller (IH-EMPC with perfect forecast) is expected to have performance bracketed between two of the previously discussed controllers. Specifically, for large horizons, it will approach the EMPC with perfect forecast, and for small horizons, it will approach the IH-EMPC with imperfect forecast. The details of IH-EMPC with perfect forecast will be presented in future work.

■

AUTHOR INFORMATION

Corresponding Author

* Tel.: 312-567-3537. Fax: 312-567-8874. E-mail: chmielewski@ iit.edu. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS D.C. would like to thank I. Grossmann and R. Lima for many insightful discussions that served as the foundation and motivation for this work. Both authors would like to thank the National Science Foundation (No. CBET-0967906) for financial support.

■

REFERENCES

(1) Robinson, P. J.; Luyben, W. L. Plantwide control of a hybrid integrated gasification combined cycle/methanol plant. Ind. Eng. Chem. Res. 2011, 50 (8), 4579−4594. (2) Jones, D.; Bhattacharyya, D.; Turton, R.; Zitney, S. E. Optimal design and integration of an air separation unit (ASU) for an integrated gasification combined cycle (IGCC) power plant with CO2 capture. Fuel Process. Technol. 2011, 92 (9), 1685−1695. (3) Mahapatra, P.; Bequette, B. W. Process design and control studies of an elevated-pressure air separations unit for IGCC power plants. In Proceedings of the American Control Conference, Baltimore, MD, June 2010; pp 2003−2008. (4) Bequette, B. W.; Mahapatra P. Model Predictive Control of Integrated Gasification Combined Cycle Power Plants: Technical Report; OSTI ID: 1026486; Department of Energy (DOE): Washington, DC, 2010. (5) Rawlings, J. B.; Amrit, R. Optimizing process economic performance using model predictive control. In Nonlinear Model Predictive Control; Lecture Notes in Control and Information Sciences, Vol. 384; Springer: Berlin, 2010; pp 119−138. (6) Rawlings, J. B. Tutorial overview of model predictive control. IEEE Control Syst. Mag. 2000, 20 (3), 38−52. (7) Braun, J. E. A comparison of chiller-priority, storage-priority, and optimal control of an ice-storage system. ASHRAE Trans. 1992, 98 (1), 893−902. (8) Morris, F. B.; Braun, J. E.; Treado, S. J. Experimental and simulated performance of optimal control of building thermal storage. ASHRAE Trans. 1994, 100 (1), 402−414. (9) Kintner-Meyer, M.; Emery, A. F. Optimal control of an HVAC system using cold storage and building thermal capacitance. Energy and Buildings 1995, 23, 19−31. (10) Henze, G. P.; Krarti, M.; Brandemuehl, M. J. Guidelines for improved performance of ice storage systems. Energy and Buildings 2003, 35, 111−127. (11) Braun, J. E. Impact of control on operating costs for cool storage systems with dynamic electric rates. ASHRAE Trans. 2007, 113 (2), 343−354. (12) Mendoza-Serrano, D. I.; Chmielewski, D. J. Infinite-horizon Economic MPC for HVAC Systems with Active Thermal Energy Storage. In Proceedings of the 51st IEEE Conference on Decision and Control; Maui, Hawaii, Dec 10−13, 2012; pp 6963−6968. (13) Karwana, M. H.; Keblisb, M. F. Operations planning with real time pricing of a primary input. Comput. Oper. Res. 2007, 34, 848−867. 3164


IGCC Power Plant Dispatch Using Infinite-Horizon Economic Model

Recommend Documents