Real Time Operational Optimization of a Seawater Desalination

Jump to Prediction of Freshwater Demand - ... optimization adjusts the feed pressure and feed flow rate according to the present freshwater demand and...
0 downloads 0 Views 1MB Size
Subscriber access provided by NEW MEXICO STATE UNIV

Process Systems Engineering

Real time operational optimization of seawater desalination system based on rolling prediction of hourly freshwater demand Quannan Zhang, Aipeng Jiang, Minliang Gong, Haokun Wang, Jian Wang, and Shu Jiangzhou Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b01449 • Publication Date (Web): 28 Jun 2018 Downloaded from http://pubs.acs.org on July 1, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 16 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Real time operational optimization of seawater desalination system based on rolling prediction of hourly freshwater demand Quannan Zhang, Aipeng Jiang*, Minliang Gong, Haokun Wang, Jian Wang, Shu Jiangzhou School of Automation, Hangzhou Dianzi University, Hangzhou, Zhejiang, 310018, China ABSTRACT: Since the optimal operation trajectory obtained under fixed freshwater demand and known operational condition couldn’t achieve cost saving when these parameters have uncertain changes along time, in this paper, a real time operational optimization method based on rolling prediction of hourly freshwater demand was proposed. First, by analyzing the historical data of the hourly water consumption in the seawater desalination system, a new method for predicting the daily water demand in the next 24 hours is proposed, and the predicted trajectory is continuously corrected and updated with present freshwater consumption. Then, with the well-established seawater desalination system model combined with the dynamics of storage tank, the real time optimal operational problem with its solving strategy was given to minimize the daily operational cost. The optimization problem with differential and algebraic equations was discretized into nonlinear programming problem by finite element collocation, and then a rolling optimization solution strategy based on simulation is used. Finally, case study was used to verify the proposed method. The results show that the proposed optimal operation method can achieve significant cost reduction, and can also overcome water level violation problem under fixed freshwater demand or conventional method. Key words: Seawater desalination; Operational optimization; Rolling prediction; Real time; Freshwater demand

1.

Introduction

Desalination is an important way to solve the shortage of clean freshwater resources. At present, the main seawater desalination methods in the world include reverse osmosis (RO), multi-stage flash (MSF) and multi-effect distillation (MED). With the advance of membrane industry and the improvement of new energy recovery devices, the reverse osmosis method has become one of the most promising desalination technologies.1-3 However, driven by high pressure, seawater reverse osmosis desalination is a high-energy consumption technique. The energy cost is the primary operational cost of the production process, and the performance of RO membrane has a great influence on water production performance. In order to ensure the performance of the reverse osmosis process and to reduce the operational cost of the full scale flowsheet, more and more attention has been paid to reduce the operational cost and investment cost through optimal operation and optimal design. To minimize the product costs, Kim et al. carefully

analyzed the mechanism and the optimization of RO process with system engineering method.4 With detailed capital and operational cost evaluation of one stage SWRO system, A hybrid optimization for design and operation was investigated by Lu et al., optimal structure and operation of the system with given seawater feed conditions were obtained. Considering the long-term effect of membrane fouling, an optimal schedule of reduce the annual water production cost, the membrane cleaning and replacement were also studied to reduce annually cost of operation and maintenance.5 Then the research is further extended to the multi-level and multi-stream system configuration design.6 Considering the variations in operating conditions (such as membrane fouling along time), Zhu et al. obtained a time-related system performance by optimizing the operation of an RO network.7 See et al. further studied the optimization of the RO network, including optimizing the RO process stages and the arrangement of membrane modules. 8 - 9 Mujtaba studied design optimization and operational optimization with different optimization

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

objectives such as maximum recovery, minimum energy consumption, and minimum pressure vessel usage. 10 Williams et al. studied the joint optimization problem of the device size and short-term operation of a single desalination system.11 Wilf and Bartels et al. considered several possible system configuration schemes and energy recovery methods to optimize the energy consumption of the RO system.12 Wilf and Klinko et al. hypothesized that the electricity price and operating conditions were constant, and then established the corresponding relationship between the minimum energy consumption and the seawater recovery rate.13 Taking into account the time-varying of seawater temperature, feed pressure and freshwater load during the operation of seawater desalination system, Choi and Kim built the full scale process model of reverse osmosis desalination system, and then analyzed the effects of operation parameters on the performance of the system.14 With the consideration of time-of-use electricity price and the buffer effect of the storage tank, Jiang and Cheng et al. obtained the optimal feed flowrate and optimal feed pressure trajectories minimizing the operational cost, computing results show that compared with conventional method, significant cost saving can be achieved by utilizing the change of electricity price and buffer effect of storage.15 Sassi and Mujtaba also considered the time variation of water load, seawater temperature and even salinity, they introduced storage tank to provide additional operational flexibility, and the optimal operation and design of the system are obtained under different conditions.16-17 Since energy consumption and energy cost of the system are not consistent with the change of the electricity price, Ghobeity and Mitsos considered the operating conditions and time-of-use price changes, and studied the optimal operation of the reverse osmosis desalination system with the lowest energy cost. The results show that this is more conducive to reduce operational cost.18 Some of the above studies have also considered the impact of electricity price on the operating costs of the system, some have taken into account the impact of storage tank and freshwater load on the operating cost of the system, and more energy saving

Page 2 of 16

potential can be realized with the consideration of the impact of electricity price and the buffer effect of the storage tank. However, the above research does not involve the prediction of freshwater demand. Since the daily freshwater demand has great volatility and uncertainty, if we can’t get a more accurate prediction of hourly freshwater demand, the optimal operation of the system can’t be realized efficiently. In this paper, takes into account the uncertainty and large range fluctuation of freshwater demand, a real time optimal operation strategy is proposed based on rolling prediction of freshwater demand. With the preliminary prediction of freshwater demand in the next 24 hours, the optimal operation takes one hour as an operational interval and is executed step by step based on freshwater demand rolling prediction and updates. Since the optimal operational trajectories are obtained by sequentially solving the dynamic optimization problem, the related computing technique is also investigated to achieve fast and stable solution.

2. Modeling and optimization of seawater reverse osmosis system The full-scale flow chart of the one stage SWRO system with high desalination membrane can be observed in Figure 1. Based on the Solution-Diffusion and Film theory, many researchers, such as Senthilmurugan and Oh, analyzed the spiral-wound membrane desalination process and established the reverse osmosis process model under different conditions. 19-23 And based on the process model simulation and optimization were conducted using the lumped parameter approximation and other computing methods. 24-25

Figure 1. Basic schematic diagram of a typical SWRO desalination system

For reverse osmosis membrane, the overall fluid and

ACS Paragon Plus Environment

Page 3 of 16 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

solute mass balance equations are as follows:

standard condition, α1, α2 and β1 are constant parameters,

Q p = Q f − Qr

(1)

Q f C f = Qr C r + Q p C p

(2)

T represents the operating temperature at Kelvin temperature. The Osmotic pressure and concentration are almost linear relationship, expressed as:

L

Q p = nlW ∫ Jvdz

∆π = RT (Cm − C p )

(3)

(10)

0

Where f, r and p refer to feed, reject (brine) and permeate (product) streams, respectively. Q and C represent the flowrate and salt concentration, respectively. The flux of water and salt in different parts can be obtained through the mass transport relationship of Muraenid - Sourirajan dissolving and diffusion. nl, W and L denote the number of leaves, the width and the length of RO module, respectively. According to the principle of solution diffusion, solvent flux and solute flux through the membrane can be expressed as: Jv = Aw (∆P − ∆π ) (4) Js = Bs (Cm − Csp )

(5)

∆P = ( Pb − Pp )

(6)

Pb = Pf − Pd

(7)

Here Aw denotes membrane water permeability, Bs denotes membrane TDS permeability, Pf denotes feed pressure, Pd is pressure of seawater in RO module channel, Pp is permeation side pressure (usually the default is the ambient pressure), ∆π denotes the osmotic pressure difference, Cm denotes solute concentration on the feed side of the membrane surface, Csp denotes solute concentration on the permeate side of the membrane surface, Cp is value of module terminal of Csp, this means that Cp=Csp (L). Aw and Bs are more sensitive with operating temperature and related factors, the relationship is as follows: T − 273 Aw = Aw0 exp(α1 − α 2 ( Pf − Pd )) (8) 273 Bs = Bs 0 exp( β1

T − 273 ) 273

(9)

where R is a gas constant. In order to solve the above equations, it is necessary to know the specific parameters of the RO process and the value of Cm. It can be derived from the steady state equilibrium and CP theory around the boundary layer.

φ=

Cm − C p Cb − C p

= exp(

Jv ) kc

(11)

The concentration of the channel Cb and the solvent transports Jv vary along the membrane channel. kc is the mass transfer coefficient, which is calculated as follows: Sh =

kc d e = 0.065Re0.875 Sc0.25 DAB

Re = ρVde / µ Sc = µ / ( ρ DAB )

(12) (13) (14)

where ρ stands for the density of permeable water, de is hydraulic diameter of feed spacer Channel, µ is kinematic viscosity and DAB is the dynamic viscosity, it can be calculated by the regression equation: DAB = (0.726 + 0.023T + 0.000277T 2 ) × 10 -9 (15)

Figure 2. Scheme of the rectangular channel model of spiral wound module feed channel

Figure 2 shows the axial variation of the feed parameters in the channel of the RO membrane. According to Figure 2, the pressure loss during RO in the membrane element is:

where Aw0 and Bs0 are the intrinsic transport parameters in

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

dPd ρ V2 = −λ dz de 2

dC p ,out

(16)

where

λ = 6.23K λ Re −0.3

(17)

Kλ denotes the empirical parameter. The bulk pressure of RO channel is Pb=Pf-Pd, so we can get: dPb dP ρ V2 =− d =λ dz dz de 2

Page 4 of 16

dt

=

Qp St H t

(Csp − C p ,out )

(24)

Here St and Ht refer to storage tank area and water level, respectively. Qout is freshwater demand for users, Ct,out is desalinated water salinity output to the user, Csp is the average salt concentration of permeate water from RO units.

(18)

Here z =0, Pb = Pf ; z =L, Pb = Pr .

H t ,up

V denotes axial velocity, and which is satisfied: dV 2 Jv =− dz hsp

Ht

(19) H t ,lo

when z =0, V = Vf =

Qf

neWhsp

Figure 3. Schematic diagram of the storage tank

(20)

when z =L, V = Vr =

Qr neWhsp

Main performance indexes of seawater reverse osmosis system include water recovery ratio, salt rejection coefficient and specific energy consumption per unit product, respectively. They can be expressed as follows:

(21)

Rec = Q p / Q f

where hsp is the height of the feed interval channel. The change of bulk concentration Cb along the membrane channel is given by the following equation: dCb 2 Jv (Cb − C p ) = (22) dz hspV And z =0, Cb = C f ; z = L, Cb = Cr . Based on the above equation, under the given manipulated variables, environmental parameters and membrane parameters, the flowrate of permeate water Qp and its concentration Cp can be obtained through solution of above equations. For the storage tank in seawater reverse osmosis system (Figure 3), negligible process of PH quenching and other treatment processes is not considered. According to mass conservation law, the dynamics of water level and salt concentration of product water is as follows: dH t = (Q p − Qout ) / St (23) dt

Ry =

C f − Cp Cf

×100% = (1 −

(25) Cp Cf

) ×100% (26)

SEC = ( Pf Q f / ε p − Pr Qr ε ef ) / Q p

(27)

Where ߝp and ߝpf are the mechanical efficiency and energy recovery efficiency, respectively. To reduce the energy cost of the system and ensure the water quality and safety requirement, the optimal operation problem can be expressed as: tf

Min

tf

∫ E p SEC Q p dt = ∫ E p ( Pf Q f / ε p − PrQr ε ef ) dt (28)

Q f , Pf , H t t 0

t0

Equality constraint: Eqn. (1) - Eqn. (27) Boundary condition constraint:

ACS Paragon Plus Environment

Page 5 of 16 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Pflo ≤ Pf ≤ Pfup

Fi,t +1 =

Q flo ≤ Q f ≤ Q fup

Fi ,t +1 = ψ 1 xi ,t + ψ 2 xi ,t -1 + ⋅ ⋅ ⋅ + ψ n xi ,t − n +1 + ψ 0

C plo ≤ C p ≤ C pup Cout ,lo ≤ Cout ≤ Cout ,up H t min ≤ H t ≤ H t max

Among them, Ep indicates the price of electricity, which is hourly different in a day. Considering the quasi- periodic changes of other parameters, such as ambient temperature, the optimal time interval is set to 0 to 24 hours. Then the objective function of the optimization problem can be expressed as: 24

Min

∑E

p

(i )SEC (i ) Q p (i )

(29)

Q f , Pf , H t i =1

Considering the consistency of water level and the safety margin, the following constraints need to be added: Ht (0) − Ht (24) ≤ ∆

(30)

Because of the buffer effect of the storage tank, the permeate water can be increased when the electricity price is low to reduce the energy cost. However, considering the safety constraints of the storage tank level and the quality restriction of the production water, the production process is limited by the water supply process based on freshwater demand. Accurate prediction freshwater demand within the next 24 hours is the precondition for optimal operation of the SWRO system.

3.

1 ( xi,t + xi ,t −1 + ⋅⋅ ⋅ + xi ,t −n+1 ) n

Prediction of freshwater demand

Without considering the storage tank and hourly electricity price change, real time operational optimization adjusts the feed pressure and feed flowrate according to present freshwater demand and present operating condition. But if we extend the optimal operation to whole day’s interval, a better prediction of freshwater demand in the future 24 hours is needed. Generally, historical water consumption can provide lots of information for our prediction. The moving average method or the autoregressive moving average method can be used for the prediction, which is as follows: 26-27

(31) (32)

Here i=1,2…24 refer to hours, xi,t represents the actual observed data at the ith hour of the tth day, n represents number of historical sample days; Ψ0…Ψn refer to the value of the coefficient obtained by regression. Figure 4 shows some historical data for an SWRO system.28 It is shown from the diagram that there is a special regularity in the freshwater load, but the fluctuation is also relatively large. Based on the moving average method or auto-regressive moving average method, it is possible to predict the next day's average water consumption and then schedule the water production process. But it is difficult to achieve the desired accuracy when the water fluctuation is large, and the cumulative effect of error on water level will cause safe problem. Based on the above considerations, this paper uses the principle of predictive control to reduce the prediction error by rolling prediction and rolling correction of above prediction. After obtaining the predicted curve for the second day, the hourly water load demand profile is updated and corrected each hour. Accordingly, with the water demand prediction, the optimal problem is solved to obtain and update the optimal operation profile (optimal feed flowrate and optimal feed pressure) each hour. And then after the optimal feed pressure and feed flowrate is executed and current water load is measured, the new procedure for water prediction and solution of the optimal problem is carried out. To get better prediction of water load, the characteristics of history data used to get the profile template, and the data from the 20th hour of the first day to the 5th hour of the second days was used to get and correct the prediction profile in a few hours ahead. The prediction procedure is as follows (as shown in Figure 5): Step 1: Obtain the hourly water load profile template with history data, the water load profile template can be calculated by moving average method as below:

Fi ,m = ( xi ,t + xi ,t −1 + ⋅⋅⋅ + xi ,t − n +1 ) / n

ACS Paragon Plus Environment

i=1,2,…,24.

Industrial & Engineering Chemistry Research

Step 2: With the hourly freshwater consumption of the first day, predict hourly freshwater demand of the first several hours in the next day by differential auto-regressive moving average method mentioned below. Then obtain the initial water load prediction profile of the whole second day. The prediction of the next hour’s freshwater load can be represented by freshwater load sequence ahead, which can be calculated by Eqn.(33). 28

Yi +1 = φ1Yi + φ2Yi -1 + ⋅⋅⋅ + φnYi − n +1 + φ0

(33)

Figure 4. Historical data of hourly freshwater consumption

Then the auto-regressive equation in differential form can be expressed as:

Yi +1 -Yi = φ1( Yi -Yi -1 ) + ⋅⋅⋅ + φn (Yi − n +1 -Yi −n ) Here

i=0, Obtain template trajectory of hourly water load Fi,m with history data

(34) Based on history data, get regression efficient through nonlinear least squares φ1 ...φn

φ1 …. φn denote the auto-regressive coefficient,

which can be obtained by nonlinear least square. With Eqn.(34) the hourly water load of the first 4-5 hour

Y25 ...Y29

is obtained sequentially.

Prediction of the first 5 hours’ water load with regression equation

And then with

comparison to profile template the whole profile prediction of the next day presented by Y25 ...Y29 ...Y48 is initially

i=i+1

Get the observed value of the ith hours’ water load of next day,and correct the predicted value

obtained. Here Y30 ...Y48 is calculated by:

Y j = F j − 24, m + Gap ,

Calculation the Gap between template trajectory and predicted load, then get the whole days’ water load prediction

j = 30...48

(35) No

j −1

Gap = Min

∑(Y

i=i+1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 16

i≥ 5

− F j − 24,m) 2

j

j = 25

(36)

Yes

th

Step 3: Correct the water load of the i hour after it is measured, and then update the hourly load from the (i+1)th hour to the 24th hour by Eqn.(35)-Eqn.(36).

Recalculate the Gap , revise the left hours’ water load prediction

No i≥ 23

Yes The end, water load prediction of the next day finished

Figure 5. Schematics of rolling prediction and correction of hourly

ACS Paragon Plus Environment

Page 7 of 16 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

K

freshwater load

4.

w( z ) = wi −1 + hi ∑ Ω q ( q =1

Whole operational optimization strategy

4.1 The whole solution framework With known water load prediction, solution of the optimal problem formulated by Eqn.(1)-Eqn.(28) leads to optimal operation profile of feed pressure and feed flowrate from the first to the 24th hour. Since the water load prediction is hourly updated, to reduce prediction errors, the optimal problem accordingly should be solved sequentially as the prediction profile changes. Since the problem in comprised with a set of differential equations and algebraic equations, it is not easy to be solved in direction way based on necessary optimality conditions. To solve this problem, simultaneous method based on finite element collocation29-31 was used to discretize the optimal problem and transfer it into nonlinear programing. Then large scale solver IPOPT was used for its efficient solution. 32 With the simultaneous method, obtaining the initial value for all discretized variables is quite critical for the quick and steady solution. So in this study, simulation based strategy was used to get reasonable value of state variables. With fixed value of control variables and initial state value on the first finite element, the simulation was carried out one by one finite element until all the state variable values obtained. Since the Eqn.(30) is not continuous derivability, it was transformed into square form as follow:

( Ht (0) − Ht (24))2 − ∆ 2 ≤ 0

(37)

The total schematic for our solution strategy is shown in Figure 6. 4.2 Problem discretization based on finite element collocation In order to solve the optimization problem of the above DAEs (differential algebraic equations), a fully discrete simultaneous solution technique based on finite element configuration is used. This method has many advantages in the Runge-Kutta discrete method,33 the state variables in the differential equation are expressed by Eqn. (38):

z − zi −1 dw ) hi d zi , q

(38)

Here wi-1 is the value of the differential variables at the beginning of element i, hi is the length of element i, dw/dzi,q denotes the value of its first derivative in element i at the collocation point q, andΩq is a polynomial of order K, satisfying Ω q (0) = 0

q = 1..., K

Ω'q ( ρ r ) = δ q,r

q = 1..., K

(39) (40)

where ߩr is the location of the ‫ݎ‬th collocation point within each element. Continuity of the differential profile is enforced by K

wi = wi −1 + hi ∑ Ω q (1) q =1

dw dzi , q

(41)

In addition, the control and algebraic profiles are approximated using a Lagrange basis representation which takes the following form: K

y ( z ) = ∑ψ q ( q =1 K

u( z ) = ∑ψ q ( q =1

z − zi −1 ) yi , q hi

(42)

z − zi −1 )ui ,q hi

(43)

Here yi,q and ui,q represent the values of the algebraic and control variables, respectively, in element i at collocation point q, and satisfying zi-1≤z≤zi, and Ψq is the Lagrange polynomial of degree K . Through the above method, the operational optimization problem of was transformed into NLP problem. By solving the NLP problem, the optimal control trajectories within 0-24 hours can be obtained. 4.3 Initial value technique based on simulation According to the overall operational optimization diagram shown in Figure 6, we need to continuously solve the discrete NLP problem. However, the solver need reasonable estimate of all variables before it start, and it is crucial to the stability and efficiency of the optimal solution. In order to obtain a reasonable initial value of all discrete variables, an initial value acquisition technique based on stepwise simulation is adopted. For the original NLP

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

problem, with given value of control variables, the original optimization problem is transformed into simulation problem with equal equality equations and variables. With initial value of state variable, simulation solution of the ith finite element can get all the values of state variable in each collocation point and thus the initial value of the (i+1)th finite element. By continually solving the simulation problem from the ith finite element to the end finite element, all the values of the discrete state variables can be obtained in reasonable way. Optimization problem formulated with equation(1-27)and(37)

Discretized to large scale NLP by finite element collocation method

Let i=0, set initial parameters; get water load prediction trajectory from (i+1)th to 24th hour

Input water load prediction trajectory from (i+1)th to 24th hour

Solve the NLP one by one finite element with given guess of control variables, until all the values of state variable obtained

Correction the water load of the ith hour, and recalculate the load trajectory from (i+1)th to 24th hour

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Free control variables and solve the NLP with large scale solver, get the optimal trajectory of control variables from (i+1)th to 24th hour

The optimal control variables at the current ith hour was manipulated for the system

At the next hour,let i=i+1

Get the value of observed water load of the ith hour

N

i≥23

Y The end, optimal operation of 0-24th hours in the next day was finished

Figure 6. Total schematics of the rolling optimal operation

5.

Page 8 of 16

Case Study

A single-stage SWRO system shown in Error! Reference source not found. was studied to demonstrate the cost saving effect with the proposed method. SW30HR-3XX membrane modules are used for salt removal. The RO unit consists of 55 parallel pressure vessels with 8 membrane modules in series in each pressure vessel. One storage tank is used to save the produced freshwater, and it works as the buffer between the water production and supply. The main structure parameters are presented in Table 1, and the hourly electricity price is presented in Figure 7. The maximum salt concentration of water production from the reverse osmosis plant is 1kg/m3, and the salt concentration of the freshwater to users is 0.8kg/m3. In order to ensure the water production to meet the freshwater demand, we adjust the water production process according to the freshwater demand. Figure 8 shows some data of historical freshwater consumption.28 It can be seen from the figure that the hourly water consumption has the obvious characteristics, there are two daily peaks of water consumption and after the second peak, and the water consumption continually drops until 5:00 am of the next day. However, due to the large fluctuation of water consumption, the prediction results based on the past few days have a large deviation from the actual values. The Figure 9 shows the deviation between the predicted values based moving average method and the actual value in the past few days. Since there are uncertainty factors, if the freshwater demand for the next 24 is just predicted once, the cumulative errors will be quite big, and will cause big problem if it used for optimal operation of RO process. Therefore, it is better for us get a quite accurate prediction profile of freshwater demand in the next 24 hours, and the optimal operation is adjusted according to the repeatedly update of the prediction profile. From the historical freshwater consumption, it is also can be found that, the first 5 hours water consumption of the next day are strongly associated with the whole day’s water consumption. Figure 10 shows a comparison of the template profile formed by historical data with the actual water

ACS Paragon Plus Environment

Page 9 of 16 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

consumption of the (n+1) day. Through the shift of the template profile, the deviation between the actual water curve on the (n+1) day and the adjusted template profile is significantly reduced. According to the freshwater demand rolling prediction method proposed in this paper, the profile of freshwater demand is renewed at different time of the next day is obtained, as shown in Figure 11. The figure only gives the prediction profile at initial time, the 5th, 12th and 18th hour and the actual water consumption. It can be seen from the figure that the relative deviation of the predicted freshwater demand at the initial time is still relatively large, but with the step-by-step update and correction, the cumulative error is getting smaller and smaller, which is as shown in Figure 12. Table 1.

Figure 7. Profile of hourly electricity price

Feed condition and tank information

Feed condition

value

No. pressure vessel

55

No. module per vessel

8 3

Feed concentration (kg/m )

30

Feed temperature (℃)

25

Feed pressure (Bar)

59

Feed PH

5-8

Tank information

value

Area of liquid level St (m2)

150

Total height (m)

17

Height warning limit (m)

15

Min height limited (m)

2

Initial concentration (kg/m3)

0.430

Initial water liquid level (m)

4

Figure 8. Hourly freshwater consumption of historical days

Figure 9. Freshwater demand prediction by average moving method

Figure 10. Comparison of real and pattern value after parallel movement

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 11. Schematic of rolling prediction of hourly water demand

Figure 12. Total error changes along time

In order to efficiently solve the operational optimization problem composed of Eqn. (1)-Eqn.(28), the differential-algebraic equations of RO module are discretized by 20 finite elements and 3 Radau collocation points; the dynamic equations of the storage tank are discretized by 24 finite elements and 3 Radau collocation points, so that the discrete precision is already relatively high. More information about the discrete optimization problem can be seen in Table 2. The discrete optimization problem was solved by large-scale NLP solver IPOPT under GAMS24.0. 34 And a 2.5GHz computer with a Win10 operating system, Intel Core i7-6500 CPU processor and 8GB of memory is used. For the purpose of energy cost saving analysis, comparison of different strategies with our proposed strategy of optimal operation was carried out. Case 1: the freshwater demand is not predicted, and the water production is adjusted according to the present water consumption and operational conditions. In this case, the freshwater consumption is consistent with the water production and only present feed pressure and feed

Page 10 of 16

flowrate are optimized. Case 2: the moving average method is used to predict the hourly freshwater demand of the second day, and once we get the optimal operational trajectories of feed pressure and feed flowrate through solving the optimization problem, the feed pressure and feed flowrate in the next 24 hours are controlled based on the above trajectories. Case 3: rolling prediction and updating of freshwater demand with our method, sequentially solving the operational optimization problem based on updated freshwater prediction, the feed pressure and feed flowrate is sequentially controlled based on present point of optimized trajectories. In order to get the reasonable initial value for the optimization solution, the simulation based technique was adopted. With our problem solving technique, total computing time of the problem is about 34.375 seconds with IPOPT. The simulation time is 23.413 seconds, and the optimal solution time is 10.962 seconds. For the successive optimization solution, the obtained optimal values of all the variables can be used as the initial value for the next computing, so the other 23 times of optimal solution of the optimization problem can be computed more stably and efficiently, the total calculation time is significantly reduced. The operational optimization problem was successfully solved under Case 1 to Case 3, calculation results are shown in Table 3 and Figure 13 to Figure 17. As can be seen from Table 3, the energy cost is quite different under the cases. In Case 1, the daily energy cost reached 9170 China Yuan. Compared with Case 1, the second case can obtain 9.49% energy cost saving, but the water level changes greatly, not only the termination level is much larger than the initial level, but also the water level is higher than the alarming level (seen in Figure 15). In the Case 3, the daily energy cost is only 6,710 China Yuan, compared with Case 1, about 26.83% energy cost saving can be achieved. Furthermore, the final water level is near the initial value, the water level constraint is satisfied. For the Case 1, the initial and final water levels are the

ACS Paragon Plus Environment

Page 11 of 16 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

same, because the hourly water consumption is consistent with the amount of water production. For the Case 2, with the total water consumption constraint ∆=0.5m, the operational problem is solved and the obtained control trajectories is used for the optimal operation. But it is found that the water level reached 13.6m, which is far from the initial value. This is cause by the accumulative error of freshwater demand prediction. For the Case 3, after successive optimal operation, the final water level, although not at the initial 7 m, satisfies the water level constraints. In Case 3, the total computing time is far greater than that of the first two Cases, this is because the operational optimization problem is solved 24 times in Case 3. Figure 13 and Figure 14 show the optimal feed pressure and the feed flowrate trajectories for the cases. Compared with the first case, feed pressure and feed flowrate of the other cases change in large quantities, because they both try to make full use of the impact of hourly electricity price change and the buffer effect of storage tank. The increase of the feed flowrate and feed pressure under the condition of low electricity caused the increase of water level, which can be seen in Figure 15. Quality of water production and water supply are two important performance indexes to be ensured, the hourly change of then under different cases can be seen in Figure 16 and Figure 17. It can be seen that, though sometimes salt concentration from RO unit is higher than 0.8kg/m3 under Case 2 and Case 3, the salt concentration for supply water is always satisfy the required demand, this is mainly due to the buffer effect of the storage tank. Table 2.

Collocation point each hour

Table 3. No. Case 1 Case 2 Case 3

3

Results of the three cases Solution time (s) 34.375 34.768 285.892

Objective (104CNY/day) 0.917 0.830 0.671

Cost save (%) 0 9.49 26.83

Final Ht (m) 7 13.6 6.68

Figure 13. Optimal feed pressure profiles of the cases

Figure 14. Optimal feed flowrate profiles of the cases

Information of the discretized problem

Model condition

value

Total Variables

31324

Equality constraints

31276

Inequality constraints

1

Nonzeros in Jacobian

10992

Nonzeros in Hessian

39626

Finite element along membrane

20

Collocation point along membrane

3

Sub-interval of time(hour)

24

Figure 15. Height of storage tank profiles of the cases

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 16. Concentration profiles of permeate water of the cases

Page 12 of 16

of related variables at different hours, which can be seen in Figure 19 to Figure 23. In these figures, profiles of variables such as feed flowrate, feed pressure at the 1, 9, 16 and 24 hours are given to show evolutionary process along time. In Figure 19 and Figure 20, the profile of the first hour indicates the optimal operational trajectories with of the original freshwater demand prediction, and the profile of the 24th hour represents the final completed optimal operation trajectories. For the optimal feed pressure profile of the 9th hour, the profile of the first 8 hours indicates that these optimal feed pressures have been executed. The subsequent optimal feed pressure represents the expected optimal feed pressure to be performed based on the 9th prediction of freshwater demand for the left 15 hours.

Figure 17. Concentration profiles of supply water of the cases

For the Case 3, due to the optimal operation is performed according to sequentially updating freshwater demand prediction and sequentially solving the optimization problem. The above results give the final optimal operational profile and the objective function value. In fact, each time of freshwater demand prediction updating and the optimization solving would result in updating of optimal operational trajectories and related state profiles. And the value of objective function changes as the updating proceeds. Figure 18 shows that the optimal energy cost obtained through optimization problem solving changes with the updating of optimal feed pressure and feed flowrate trajectories. Figure 19 and Figure 20 show the continuously updating of optimal feed pressure and feed flowrate curves. Since the freshwater demand prediction and correction is executed in each hour until the next 24th hour, through 24 times of optimization problem solving and control trajectories tracing, we can get 24 different profiles of control variables as well as state variables and related performance index. For simplicity, we just give the profiles

Figure 18. Changes of optimal energy cost with hourly updating

Figure 19. Rolling update profile of optimal feed flowrate

ACS Paragon Plus Environment

Page 13 of 16 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 20. Rolling update profile of optimal feed pressure

Figure 21. Rolling update profile of storage height

Figure 22. Rolling update profile of supply water concentration

Figure 23. Rolling update profile of permeate water concentration

6.

It is of great significance to reduce energy cost of SWRO system by adjusting the feed pressure and the feed flowrate according to the change of operating conditions. But if we just consider the water production process and just minimize specific energy consumption, the energy cost saving potential is limited. Considering the freshwater demand of users and water production process, utilizing the storage tank and hourly electricity price to add operational flexibility can significantly improve the energy cost potential, but this requires accurate freshwater demand information and suitable operation strategy. In this paper, with the consideration of linkage between water production and water supply process, a real time optimal operation strategy based on rolling prediction of freshwater demand is proposed to reduce daily energy cost. Based on historical data analysis and the idea of predictive control, a better freshwater demand rolling prediction is obtained, which is used to predict and rolling updating the freshwater demand profile of the next 24 hours. With the predicted freshwater demand profile, the optimal operation problem is sequentially solved, leading to the new optimal operational trajectories. By sequentially control the feed pressure and feed flowrate to the present point of the optimal operational trajectories, the real time optimal operation is realized. For efficient solving the operational optimization problem and realization of the proposed method, some computing techniques are also studied. Computing results of a typical seawater desalination system under different optimal operation cases show that: Compared with the non-scheduling operational optimization mode, the daily energy cost of the system can be significantly reduced by 26.83% with the proposed strategy. And compared with the operational optimization mode based on freshwater demand prediction of moving average method, the energy cost with the proposed strategy can also reduce significantly, and at the same water level constraint for safety can also be well guaranteed.

Author Information

Conclusions

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Corresponding Author *Email: [email protected]. Tel:±86 13968138854. Fax:±86 571 86919133. Notes The authors declare no competing financial interest.

Acknowledgments This work was supported by the National Natural Science Foundation of China (No. 61374142), the Public Projects of Zhejiang Province, China (No. 2017C31065) and the Natural Science Foundation of Zhejiang (LY16F03006).

References [1] Cameron, I.; Clemente, R. SWRO with ERI's PX pressure exchanger device–A global survey. Desalination, 2008, 221 (1-3):136-142. [2] Greenlee, L.; Lawler, D.; Freeman, B.; Marrot, B.; Moulin, P. Reverse osmosis desalination: water sources, technology, and today's challenges. Water Research, 2009, 43(9):2317-2348. [3] Ghobeity, A.; Mitsos, A. Optimal design and operation of desalination systems: new challenges and recent advances. Current Opinion in Chemical Engineering, 2014, 6:61-68. [4] Kim, Y. M.; Kim, S. J.; Kim, Y. S.; Lee, S.; Kim, I. S.; Kim, J. H. Overview of systems engineering approaches for a large-scale seawater desalination plant with a reverse osmosis network. Desalination, 2009, 238 (1-3):312-332. [5] Lu, Y. Y.; Hu, Y. D.; Xu, D. M.; Wu, L. Y. Optimum design of reverse osmosis seawater desalination system considering membrane cleaning and replacing, Journal of Membrane Science, 2006, 282 (1-2):7-13. [6] Lu, Y. Y.; Hu, Y. D.; Zhang, X. L.; Wu, L. Y.; Liu, Q. Z. Optimum design of reverse osmosis system under different feed concentration and product specification. Journal of Membrane Science, 2007, 287(2): 219-229. [7] Zhu, M.; El-Halwagi, M.; Al-Ahmad, M. Optimal

Page 14 of 16

design and scheduling of flexible reverse osmosis networks. Journal of Membrane Science, 1997, 129(2): 161-174. [8] See, H.; Vassiliadis, V.; Wilson, D. Optimization of membrane regeneration scheduling in reverse osmosis networks for seawater desalination. Desalination, 1999, 125 (1-3):37-54. [9] See, H.; Wilson, D.; Vassiliadis, V.; Parks, G. Design of reverse osmosis (RO) water treatment networks subject to fouling. Water Science and Technology, 2004, 49(2): 263-270. [10] Villafafila, A.; Mujtaba, I. M. Fresh water by reverse osmosis based desalination: Simulation and optimization. Desalination, 2003, 155 (1): 1-13. [11] Williams, C. M.; Ghobeity, A.; Pak, A. J.; Mitsos, A. Simultaneous optimization of size and short-term operation for an RO plant. Desalination, 2012, 301:42-52. [12] Wilf, M.; Bartels, C. Optimization of seawater RO systems design. Desalination, 2005, 173(1):1-12. [13] Wilf, M.; Klinko, K. Optimization of seawater RO systems design. Desalination, 2001, 138(1-3): 299-306. [14] Choi, J. S.; Kim, J. T. Modeling of full-scale reverse osmosis desalination system: Influence of operational parameters. Journal of Industrial and Engineering Chemistry, 2015, 21 (2):261-268. [15] Jiang, A.; Cheng, W. Operational optimizations of full flowsheet spiral-wound seawater reverse osmosis system. Journal of Chemical Industry, 2014, 65(4):1133-1143. [16] Sassi, K. M.; Mujtaba, I. M. Effective Design of reverse osmosis based desalination process considering wide range of salinity and seawater temperature. Desalination, 2012, 306 (18) :8-16 [17] Sassi, K. M.; Mujtaba, I. M. Optimal operation of RO system with daily variation of freshwater demand and seawater temperature. Computers & Chemical Engineering, 2013, 59 (4): 101-110 [18] Ghobeity, A.; Mitsos, A. Optimal time- dependent

ACS Paragon Plus Environment

Page 15 of 16 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

operation of seawater reverse osmosis. Desalination, 2010, 263(1–3): 76-88 [19] Senthilmurugan, S.; Ahluwalia, A.; Gupta, S. K. Modeling of a Spiral-wound module and estimation of model parameters using numerical techniques. Desalination, 2005, 173(3): 269-286 [20] Marriott, J.; Sorensen, E. A general approach to modeling membrane modules. Chemical Engineering Science, 2003, 58(22): 4975-4990 [21] Kaghazchi, T.; Mehri, M.; Ravanchi, M. T.; Kargari, A. A mathematical modeling of two industrial seawater desalination plants in the Persian Gulf region. Desalination, 2010, 252 (1-3): 135-142 [22] Li, M. Optimal plant operation of brackish water reverse osmosis (BWRO) desalination. Desalination, 2012, 293:61-68 [23] Wijmans, J. G.; Baker, R. W. The solution- diffusion model: a review. Journal of Membrane Science, 1995, 107(1-2): 1-21 [ 24 ] Boudinar, M. B.; Hanbury, W. T.; Avlonitis, S. Numerical simulation and optimization of spiral-wound modules, Desalination, 1992, 86(3): 273-290 [ 25 ]Sassi, K. M.; Mujtaba, I. M. Simulation and optimization of full scale reverse osmosis desalination plant. Computer Aided Chemical Engineering, 2010, 28: 895-900 [ 26 ]Zhang, X.; Dang, Z.; Zhang, X.; Ding, M. Comprehensive study of urban water consumption prediction model. Journal of Water Resources and Water Engineering, 2005, (04): 28-32 [27]Wu, Shao-chi.; Wu, Y.; Xu, P. An improved time series segmentation trend forecast method. Journal of Electronic Measurement and Instrumentation, 2010, (10): 953-957. [28] Jiang, Z. Research and development on optimal operation and scheduling system in desalination. Hangzhou: Hangzhou Dianzi University, 2012 [29] Biegler, L. T. An overview of simultaneous strategies for dynamic optimization. Chemical Engineering and

Processing, 2007, (46):1043–1053 [30] Bartl, M.; Li, P.; Biegler, L. T. Improvement of state profile accuracy in nonlinear dynamic optimization with the quasi-sequential approach. AICHE Journal, 2010, 57(8): 2185-2197. [31]Jiang, A.; Wang, J.; Cheng, W.; Xing, C.; Jiangzhou, S. A dynamic optimization strategy for the operation of large scale seawater reverses osmosis system. Mathematical Problems in Engineering, 2014, (2014):1-12 [32] Biegler, L. T. Nonlinear Programming: Concepts, Algorithms and Applications to Chemical Processes. Pittsburgh: Society for Industrial and Applied Mathematics (Cambridge University Press), 2010 [33] Biegler, L. T.; Zavala, V. M. Large-scale nonlinear programming using IPOPT: An integrating framework for enterprise-wide dynamic optimization. Computers & Chemical Engineering, 2009, 33(3): 575-582 [34] GAMS Development Corporation, https//:www.gams.com

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Optimal flowrate profile of seawater was obtained with predicted hourly freshwater demand profile, and it was sequentially updated according to updated prediction of hourly freshwater demand.

For Table of Contents Only

ACS Paragon Plus Environment

Page 16 of 16