IN MEMBRANE DISTILLATION PLANTS

CROSS REFERENCE TO RELATED APPLICATIONS

[0001] This application claims priority to, and the benefit of, co-pending U.S. provisional application entitled "NONLINEAR LYAPUNOV-BASED BOUNDARY CONTROL OF

DISTRIBUTED HEAT TRANSFER MECHANISMS IN MEMBRANE DISTILLATION PLANTS" having serial no. 62/187,618, filed July 1 , 2015, which is hereby incorporated by reference in its entirety.

BACKGROUND

[0002] A membrane distillation (MD) process is thermal, and includes a transfer of water vapor from a feed container to a permeate container through a micro-porous membrane, and reject other non-volatile constituents present in the influent water. The process is driven by the temperature difference along the membrane boundary layers. This gives significance to the temperature distribution among the process domain and especially along the membrane boundary layers. However, the MD process has a serious drawback in the production rate stability, and its energy consumption.

SUIVIIVIARY

[0001] Embodiments of the present disclosure are related to control of membrane distillation (MD) processes.

[0002] in one embodiment, among others, a system comprises a membrane distillation (MD) process comprising a feed side and a permeate side separated by a membrane boundary layer; and processing circuitry configured to control a water production rate of the D process based at least in part upon a distributed heat transfer across the membrane boundary layer, in one or more aspects of these embodiments, the processing circuitry can maintain a temperature difference along the membrane boundary layer at a level that promotes the water production rate of the D process. Control of the water production rate can be based upon nonlinear Lyapunov-based boundary control of the temperature difference of membrane boundary layer, perturbation-based extremum seeking control (PESC) of the temperature difference of membrane boundary layer, or Newton-based multivariabie extremum seeking control (ESC) of the temperature difference of membrane boundary layer. In one or more aspects of these embodiments, a semi-discretized model of the membrane boundary layer can be utilized. The processing circuitry can monitor inlet temperatures of the feed side and the permeate side to control the water production rate of the MD process. The inlet temperatures of the feed side and the permeate side can be bounded by a defined minimum temperature and a defined maximum temperature. The processing circuitry can monitor inlet flow rates of the feed side and the permeate side to control the water production rate of the D process.

[0003] In another embodiment, a method comprises determining a plurality of estimated temperature states of a membrane boundary layer separating a feed side and a permeate side of a membrane distillation ( D) process; and adjusting inlet flow rate or inlet

temperature of at least one of the feed side or the permeate side to maintain a difference temperature along the membrane boundary layer about a defined reference temperature based at least in part upon the plurality of estimated temperature states. In one or more aspects of these embodiments, control of a water production rate of the D process can be based upon nonlinear Lyapunov-based boundary control of the difference temperature of membrane boundary layer, perturbation-based extremum seeking control (PESC) of the difference temperature of membrane boundary layer, or Newton-based multivariabie extremum seeking control (ESC) of the difference temperature of membrane boundary layer. The difference temperature along the membrane boundary layer can be determined based upon temperature estimates generated by a non-linear observer. The method can comprise adjusting a combination of inlet flow rate and inlet temperature of the feed side or the permeate side to control the water production rate of the D process. [0004] Other systems, methods, features, and advantages of the present disclosure will be or become apparent to one with skill in the art upon examination of the following drawings and detailed description, it is intended that all such additional systems, methods, features, and advantages be included within this description, be within the scope of the present disclosure, and be protected by the accompanying claims. In addition, all optional and preferred features and modifications of the described embodiments are usable in all aspects of the disclosure taught herein. Furthermore, the individual features of the dependent claims, as well as all optional and preferred features and modifications of the described embodiments are combinabie and interchangeable with one another.

BRIEF DESCRIPTION OF THE DRAWINGS

[0005] Many aspects of the present disclosure can be better understood with reference to the following drawings. The components in the drawings are not necessarily to scale, emphasis instead being placed upon clearly illustrating the principles of the present disclosure. Moreover, in the drawings, like reference numerals designate corresponding parts throughout the several views.

[0006] FIG. 1 is a schematic diagram depicting an example of a MD process in accordance with various embodiments of the present disclosure.

[0007] FIG. 2 is an example of discretization of a feed subsystem and a permeate subsystem of the MD process of FIG. 1 in accordance with various embodiments of the present disclosure.

[0008] FIG. 3 is an example of a schematic diagram for the grid of discretization with boundary conditions including boundary conditions in accordance with various embodiments of the present disclosure.

[0009] FIGS. 4A and 4B illustrate the convergence and error analysis of a nonlinear Lyapunov-based observer in accordance with various embodiments of the present disclosure. [0010] FIG. 5 is a block diagram illustrating an example of the structure of a nonlinear observer-based Lyapunov controller in accordance with various embodiments of the present disclosure.

[0011] FIG. 6 is a table illustrating an example of MD process parameters in

accordance with various embodiments of the present disclosure.

[0012] FIG. 7 is a plot illustrating the boundary layer temperature difference along the membrane module height in an open loop MD process in accordance with various embodiments of the present disclosure.

[0013] FIG. 8 includes surface plots illustrating examples of the behavior of the MD process at steady-state among the process domain in accordance with various embodiments of the present disclosure.

[0014] FIGS. 9A and 9B include plots illustrating examples of the temperature difference in different locations along the membrane module length in accordance with various embodiments of the present disclosure.

[0015] FIGS. 0A through 10C are plots illustrating the evolution of control law inputs and the evolution of the reference tracking error of the MD process with time in accordance with various embodiments of the present disclosure.

[0016] FIG. 10D includes surface plots illustrating examples of the behavior of the MD process at steady-state with the Lyapunov control in accordance with various embodiments of the present disclosure.

[0017] FIGS. 11 A and 11 B include plots illustrating the effect of the controller gain on the control law in accordance with various embodiments of the present disclosure.

[0018] FIG. 12 is a schematic diagram illustrating an experimental set-up in accordance with various embodiments of the present disclosure.

[0019] FIGS. 13A, 13B and 15 are tables indicating characteristics of the membrane and feed water used in the experimental set-up of FIG. 12 in accordance with various embodiments of the present disclosure. [0020] FIGS. 14A-14C and 16A-16D illustrated experimental and modeled results of the D process of FIG. 12 in accordance with various embodiments of the present disclosure.

[0021] FIGS. 17A through 17C are schematic diagrams illustrating examples of perturbation-based extremum seeking control (PESC) in accordance with various embodiments of the present disclosure.

[0022] FIGS. 18A and 18B illustrate the convergence and error of the feed and permeate states with PESC in accordance with various embodiments of the present disclosure.

[0023] FIG. 19 is a block diagram illustrating an example of the structure of a PESC controller in accordance with various embodiments of the present disclosure.

[0024] FIG. 20 includes a table indicating characteristics of the membrane used in simulations of the PESC in accordance with various embodiments of the present disclosure.

[0025] FIGS. 21A through 21C illustrates simulation results of PESC implementations in accordance with various embodiments of the present disclosure.

[0026] FIG. 22 is a schematic diagram illustrating an example of a MD cell in accordance with various embodiments of the present disclosure.

[0027] FIG. 23 is a schematic diagram illustrating an electrical analog circuit for a MD module including the MD ceils of FIG. 22 in accordance with various embodiments of the present disclosure.

[0028] FIG. 24 is a schematic diagram illustrating an example of a solar powered MD process in accordance with various embodiments of the present disclosure.

[0029] FIGS. 25A and 25B illustrate simulation results for the MD process of FIG. 24 in accordance with various embodiments of the present disclosure.

[0030] FIG. 26 is a schematic diagram illustrating an example of Newton-based multivariable ESC in accordance with various embodiments of the present disclosure. [0031] FIG. 27 is a table indicating characteristics of a D module used in the simulation of the MD process of FIG. 24 in accordance with various embodiments of the present disclosure.

[0032] FIGS. 28A and 28B are simulation results for the Newton-based multivariable ESC in accordance with various embodiments of the present disclosure.

[0033] FIG. 29 is a schematic block diagram that illustrates an example of processing circuitry employed in the membrane distillation (MD) process of FIG. 1 in accordance with various embodiments of the present disclosure.

DETAILED DESCRIPTION

[0034] Disclosed herein are various embodiments related to boundary control in membrane distillation (MD) processes. Four common configurations include air-gap membrane distillation (AGMD), vacuum membrane distillation (VMD), sweeping-gas membrane distillation (SGMD), and direct-contact membrane distillation (DCMD). Control of MD processes can help with potential commercialization of MD systems, and take advantage of their interesting properties. For example, nonlinear Lyapunov-based boundary control can be applied to distributed heat transfer mechanisms in membrane distillation plants. Reference will now be made in detail to the description of the embodiments as illustrated in the drawings, wherein like reference numbers indicate like parts throughout the several views.

[0035] Boundary controls for the temperature difference of membrane distillation boundary layers are presented. The heat transfer mechanisms inside the MD process (e.g., DCMD) are modeled with a 2D advection diffusion equation (ADE). The model can be semi- discretized in space, and a nonlinear state-space representation can be provided. The control can be designed to force the temperature difference along the membrane sides to track a desired reference asymptotically, and hence a desired flux would be generated. Certain constraints can be put on the control law inputs to be within an economic range of energy supplies. The effect of the controller gain will be discussed. Simulations with real process parameters for the model, and the controller are provided.

[0036] To achieve the desired production rate for an IV3D plant, the temperature difference along the membrane boundary layer can be controlled, because of its driven role in the MD plant. The control strategies include considering both the feed and permeate inlet temperatures as control inputs leading to a boundary control. The heat transfer mechanisms that take place in the membrane distillation process can be modeled by a 2D advection diffusion equation. A Lyapunov-based boundary control for a nonlinear process, such as the MD, is presented. Controlling partial differential equations (PDE) models can be applied over the continuous model, a semi-discretized model, or a fuily-discretized model. In this disclosure, a semi-discretized model is considered. The continuous PDE model is transformed from an infinite dimension system into a high dimension system of nonlinear ordinary differential equations (ODEs) given in the state-space representation. The advantages of adopting the semi-discretized model include that a good approximated model can be granted by refining the grid, and tuning the spatial step sizes. Furthermore, standard control strategies applicable for state-space representations can be extended to such systems.

PROCESS MODELING

[0037] The dynamics of the heat transfer accompanied in the process of membrane distillation (MD) is modeled using the 2D advection diffusion equation (ADE). The convection mechanism can be assumed to be along the z-axis (membrane length), where the diffusion mechanism can be assumed to be along the x-axis (normal to the membrane), A further extension can be applied over the model, which can lead to better consistency and flexibility in manipulating the equation by an appropriate choice of the values of the flow rates and diffusion constants, in addition, the heat convection and heat diffusion can be extended along the x and z-coordinates, even if it is for a small amount. [0038] Referring to FIG. 1 , shown is a schematic diagram depicting an example of a MD process 100. The counter current MD process 100 includes a feed subsystem 103, a membrane 106, and a permeate subsystem 109. The MD process model considers the module of the membrane to be counter current for the hot and cold water streams. The heat transfer between the feed and the permeate subsystems 103 and 09 is driven by the transmembrane heat flux, and the latent heat of evaporation. Both of them are functions of the temperature difference between the membrane boundary layers of the feed and the permeate subsystems 03 and 109. This temperature difference makes the model nonlinear.

[0039] Mathematical model. The extended model of the MD process is general.

Moreover, it is symmetric for the convection and the diffusion heat transfer mechanisms, which take place in the process domain. The model considers computing the time evolution of the temperature. In the system domain, temperature T(x, z, t) is the variable of interest with feed inlet and permeate inlet as process inputs. The mathematical model for the heat transfer mechanisms in the MD process is illustrated in FIG. 1. A rectangular domain has final values for the time and the space coordinates as shown below. The model is based on the following 2D advection diffusion equation.

,z, t) d ² T(xgj) ⁽¹⁾ with 0 < t < Tfinai , 0 < x < X, 0 < z < Z, where x e [o, x _m ] is in the feed and x e [x _pm , A']d is in the permeate. The terms v _x and v _z are the hot and the cold water stream flow rates along the x and z axes, respectively, which are assumed to be constant along the operation of the process. The terms a _x and a _z are constants that depend on physical properties of the used water: =—, (2) pep where κ _β is the thermal conductivity, c _p is the specific heat, and p is the density. The initial conditions can be set to be constant throughout the model domain as given by:

Tj{x,z, 0)™ Tf _initial f° ^{r¾ e} [°>*f _m ]» T-p(x,z,Qi) = Tp for x e [x _pm ,x] (3)

[0040] The associated boundary conditions with the model (e.g., feed solution) are assumed to be a Dirichlet condition for the inlet of each subsystem container as given by: T(x,Q,t) = T _0f forx e [ ,x _fm ], ₍₄₎

T(x,Z,t) = T _Qp forxe{x _pm ,X] _{y (5)} and Neumann conditions for the rest of the boundaries, such as: the process is isolated from the left side boundary of the feed side as given in: dT(x,z ) _!

0. (6) dx

and the right boundary side for the permeate side as given in

^' (x.z.t) and has a trans-membrane heat transfer formula for the heat flux and the latent heat of evaporation with the other subsystem (e.g. the permeate subsystem 109) such as given in:

¾ ^ _A . = [JiT)H(T) + ^k nx _lm ,zJ)

(8) for the feed side and as given in: dT(x _tZ ,t),

- {J(T)H(T) + T(x _fm , _Z: .t

(9)

for the permeate side, and finally the process has a heat flux release at the bottom boundary for the feed side as given in: dz,

and at the top boundary for the permeate side as given in:

The term k _m is the average thermal conductivity of membrane, the term S _m is the membrane thickness, φ e M ⁺ is a positive constant that represents the amount of flux that is released outside the process domain (or containers),. It may be considered to be zero in the simulations. /(T) is the trans-membrane heat flux, and H (T) is the latent heat of evaporation, which is responsible for the released sensible heat. The empirical formula for the latent heat of water is: i-J(T') = 2500,8 - 236Γ + 0.0016'f ² - 0.00006Γ 9

The trans-membrane heat flux can be based on a Knudsen diffusion model for the mass transfer coefficient given by:

where r is the membrane pore size, M is the molecular weight, R is the gas universal constant, χ is the membrane tortuosity, and e is the membrane porosity. The term AP(T) is the vapor pressure gradient between the membrane sides, which is a function of the vapor temperature based on the Antoine equation. Both are nonlinear functions of the temperature distribution.

[0041] Semi-diserefized model. The space discretization can be done using finite difference discretization schemes. Based on this, the PDE can be transformed from infinite dimensional problem into a finite dimensional group of ODEs. The convective and diffusive terms can be solved with a first-order upwind scheme and a central differencing scheme, respectively. The upwind scheme warrants the stability and ensures the solution to satisfy the entropy condition. The central difference scheme is suitable for second derivative diffusive terms. Referring to FIG. 2, shown is an example of the preparation of discretization for the feed subsystem 03 and the permeate subsystem 09 of FIG. 1 . Each point in the grid 203 depends on itself and the surrounding four points as shown in FIG. 2. Distances between two consecutive points in the same row in the grid is equal to the horizontal axis step size Ax, where it is equal to Az between two points below each other in the vertical axis. A schematic diagram for the grid of discretization with boundary conditions is shown in FIG. 3. The discretization starts in a row-wise fashion solving for the temperature of the feed subsystem 103 (FIG. 1), then solving for the temperature of permeate subsystem 109 (FIG. 1). The difference between the boundary layer ceils drives the trans-membrane heat flux and runs the D process.

[0042] The discrete version of the continuous 2D advection diffusion equation (ADE) is a nonlinear high dimension system of ODEs. The discrete model is shown in:

T(x,z,t) = AT(x, _Zi t) +BU(t) +b{T{x _t z,t)) _f (12) where vector T e Μ ^{2η χ1} is the state vector which represents the temperature distribution (at each grid point) inside the feed and the permeate domain (Tf, Tp) for the time evolution, such as in:

T T-f thus Tf, Tp 6 ^{n ' X i} with n the number of discretization points in x and z. Matrix A e I represents the operator matrix that has the dynamics of the feed and the permeate subsystems 103 and 109 with their couplings on the boundaries as shown by:

A where A _f e R" ² *" ² , , Cf _p E C _Pf e K' ¹ ' *' ¹ ' are given in:

(17)

Ax ² S _m k(

The terms and k _p are the thermal conductivity constant of seawater and the fresh water, respectively, the term 8 _m is the membrane thickness, and term k _m is the average thermal conductivity of membrane and vapor.

[0043] Matrix B e Μ ^2η2χ2 is the input matrix as shown in:

U(t) is the input vector that contains feed and permeate temperature inlets as given in:

Vector h(T) e IE ² " ^xl contains the boundary conditions associated with the membrane boundary layer cells for the feed and the permeate subsystems 103 and 109, such that the trans-membrane heat flux and the latent heat of evaporation. Both of them are functions of the membrane boundary layers temperature difference, such as in:

Vectors ζ and η in

represent the polynomial nonlinearity structures in the boundary vector. The measured output can be given by y = CT with an output matrix of C = [0 ··· 0 1 1 0 ··· 0] e M ^lx2n

[0044] Nonlinearity Analysis. The MD dynamic model is nonlinear, and is originated from the functions J( ) and H(T). However, the present nonlinearity is Lipshitz continuous, and can be proved by examining each function separately. Consequently, H(T is Lipshitz continuous following the fact that it is a smooth and infinitely differentiable polynomial function. In the other hand, !(T) can be defined as:

where c _1; c ₂ are constants, and T e [a, b], b > a > 0. Therefore, J(T is continuous and differentiable on this interval [a, b] , and bounded < M. As a result, there exists a

Lipshitz constant λ that satisfies || (Τ _Χ ) - /(T ₂ ) || < || Ti - T ₂ || and the nonlinearity in the MD dynamic model is Lipshitz continuous. CONTROLLER DESIGN

[0045] Nonlinear systems are challenging in terms of stability and control, due to their highly individualistic nature, and high sensitivity to initial and forcing conditions. The operation of heat transfer, inside the MD module, utilizes a sufficient temperature difference along the membrane boundary layers length. This boundary temperature difference is the driving force for the operation, and it causes the heat flux to be transferred from the feed subsystem 103 to the permeate subsystem 09 across the membrane 106 of FIG. 1 .

[0046] The MD process could be efficient, if the production rate of the clean water is maintained almost constant with relative low energy consumption. This can be done if a controller keeps the temperature difference of the membrane boundary layer within an effective specific range, and the inputs within economic energy supplies.

[0047] A Lyapunov-based boundary controller can be used, in order to control the boundary layers temperature differences through the inlet temperatures (and/or the inlet flow rates) for each feed and permeate subsystems 103 and 09. Satisfying the stability conditions of the Lyapunov function guarantees boundedness of trajectories, and hereby the asymptotic convergence of the states. The control law can be limited within a sufficient temperature range for heat flux production, at the same time an economic range is limited for the energy consumption. For such criteria the temperature of both feed and permeate inlets should not fall below 20°C (t^un) and not to exceed 60°C (U _j nax)-. ^as shown in: < U _min < Wf) < Umax-, Vi > 0. ₍ 23) The M D model can be rewritten to be affine system as given in:

Hence f(T) is: f(Z) =AT + b(T). (25) The output of the model y e R ^nxl is a linear combination of the difference for the boundary layer temperatures through the output matrix C, where C E m. ^nx2nI as given in: 0 ^■ · 1 0 _1 .. . 0 ■· · 0

0 · ^{■ ■} - 0 1 . .. - 1 0 . .. o

C (26)

0 ^■ · ^■ · 0 0 ^" · . 0 0 0

0 ^{■■ ■} · 0 0 1 -1 0

[0048] Proposition 1: Assume that the temperature difference of the boundary laye obey equations (24) and (25). If the following control input is used:

U = (CB) ^{~ 1} {-CAT - Cb(T) - kiCT - j). where k e JR ⁺ , then the closed loop system exhibits asymptotic reference tracking; lim where the difference between the output of the model y and the desired reference >y is defined as the reference tracking error, such as given in: e{t) = y(t) - y _r (t) . ₍₂₉₎

[0049] When the reference tracking error is asymptotic stable, the convergence of the states to their equilibrium set can be guaranteed, and so the output should track the desired reference. The analysis of the stability for the reference tracking error function is utilized by the following Lyapunov function:

V (30) (t) ^T e{t) , (31) for which the derivative is expressed by: V(t) = (t f

(CT - v (CT - ) e(i) ^T C(f(t) + BU(t)) (32) in the design stage of the Lyapunov-based boundary controller, the derivative of the Lyapunov function in equation (32) should be negative definite for global stability. The derivative of the Lyapunov function is imposed to be equal a negative value, such as given in: V (t) = -ke , (33) where k ε !E ⁺ is a positive scalar controller gain. Hence,

e(t) ^T (CAT + Cfc(r) + CB) = -ke(t) ^T e(t)■ (34) After, an adequate control law input U_(t) is selected to satisfy (34) such that the controi law is defined as:

U = - k{CL-y VJ)) _t (35) where applying this control law in equation (35) to equation (34) achieves the asymptotic stability for the predefined reference tracking error, and track the desired reference. Hereby, the derivative of the Lyapunov function is negative definite.

[0050] State transformation. Consider the following change of variables:

T = 9£Ί , (36) where T £ M ^nxl is the temperature difference vector along the membrane boundaries.

Ή E Wi $nx2n ² s a non-square transformation matrix that forms the difference of each corresponding boundary, such as

0 1 0 - 1 0 0

0 0 i 1 0

(37)

0 0 0 0

0 0 ■■ 1 0

To overcome any possible numerical issues, consider the right pseudo inverse M ^'"1 = Ή ^Ύ■ (W ^■ H ⁷ ) ^{" 1} , which as the inverse matrix of the non-square transformation matrix K, and this is because if has low condition number. By substituting equation (36) in equation (12), the structure for the MD dynamic model can be given as: where Λ = HAM ^"1 e M ^nxn represents the dynamics of the membrane boundaries.

Έ = 0ΪΒ £ R ^nx2 is the input matrix. Τ(Ή: ^{~ 1!} Τ) = Kf(T) e R ^nxl is the nonlinear effect on the membrane boundaries where f(T) = b of equation (20), and C = CM ^""1 e R ^lxn is the output matrix as defined earlier. Accordingly, the structure of the D process has been reduced, and therefore, will be easier to examine and guarantee the controllability and the

observability of the states on the membrane boundaries through the pairs {Λ;Έ) and (<A, C), respectively. Numerical simulations have shown that the controllability and the observability conditions for the MD process are satisfied.

Following Proposition 1 above, if the following control law is applied:

U - ^~% (~*.(f ~ y _r } ~ AT ~- FiV^T) , (39) where κ e M ⁺ is the controller gain, then:

e{t (AT + U{t) + {K- ^l f)) -→e(tfe(t) . (40) After, the control law will be defined as equation (39), hence asymptotic stability for the error and reference tracking is guaranteed.

NONLINEAR OBSERVER

[0051] A nonlinear observer can be used faking advantage of the well-modeled nonlinearity structure in the process. However, it can be extended to include an online constructive algorithm for finding an adequate observer gain. Consider the nonlinear Lyapunov-based observer in the following form: f - AT - BU - Τ(Ή. ^"'ι } + L(y - CT) , (41) where T 6 Ε ^{2η2 χ1} is the estimated state vector, and L is the observer gain matrix.

Following, the observer error dynamic equation can be formulated as:

T— (A - LC) - (n ^~l t) ÷ i ^ ) , (42) where ™ T ·- T is the state estimation error. The following provides sufficient conditions for exponential stability of the state observer error of equation (42), [0052] Proposition 2: There exists an exponential stable observer of the form of equation (41), if and only if, there exists a symmetric positive definite matrix, P, such that the following standard algebraic Ricatti equation is satisfied:

ΛΡ 4- PA ^r■■ ! ^■■ PRiP 4 Q, - 0 _{; (43)} where

Ri /(« 4- 1), e > 0, e€ E ⁺ - i ⁴⁴ ) and the observer gain matrix L = ^ ³ /^-

[0053] The performance of the nonlinear observer can be tested with the convergence time and error. FIG. 4A illustrates the fast convergence of the estimated states using the nonlinear Lyapunov-based observer for different locations along the membrane length. These locations are marked as a percentage of the total length of the membrane, starting from the top of the membrane to the bottom, therefore the plot shows the observer convergence on the top (6%), the middle (34%) and the bottom (67%) of the whole membrane length. The fast convergence of the observer is guaranteed by taking a smaller sampling time, in the same context, the error analysis depicted in FIG. 4B illustrates the exponential convergence of the error to zero.

NONLINEAR OBSERVER BASED LYAPUNOV CONTROLLER

[0054] Following the design of the nonlinear Lyapunov-based observer, the temperature difference along the membrane boundaries is available at each instant of time. Proceeding further, the estimated states can be feed back to a nonlinear Lyapunov controller to operate, which is valid because the observer convergence is fast. Combining the control law in equation (39) with the estimated temperature difference along the membrane boundaries in equation (41), will give the control law:

U(t) = Β ^~ "(-/ Τ - j _r ) - ΛΎ - TiH^f).- (45) FIG. 5 illustrates an example of the structure of the nonlinear observer-based Lyapunov controller.

[0055] The convergence of the combined system can be examined by studying the stability of the augmented system of the dynamics of the tracking reference error e, and the observer error T. The augmented error matrix is shown in:

and can be further simplified to be: ι? = Α¾Η (¾ , Ϊ, Τ), (47) where

represent the error state vector and the dynamics of the error, respectively. The nonlinear terms are represented in:

Sufficient conditions for asymptotic stability for the augmented error matrix are given by the following proposition.

[0056] Proposition 3: The augmented error matrix in equation (47) is asymptotic stable, if and only if, there exists matrix S = S ^T > 0, and a positive scalar ε such that the following linear matrix inequality (LMI) is satisfied:

The following shows the asymptotic stability for the error dynamics of the combined system (T→ 0, e→ 0). It can be shown that if A' and Y are real vectors of the same dimension, then for any scalar ε > 0, the following inequality holds: X ^T Y ^■■ [ ^■■ F ^r < £ ^Γ 4- (49) Consider the following Lyapunov function candidate:

V - ΰ ^τ $ύ. S - ^! > 0■ (50) Substituting the time derivative of V along the solution trajectories of equation (47) into equation (50) obtains:

= [At? 4- ®{H ^" f _} Ύ)] ^τ 8-ΰ 4- i> ^T S{&& 4- Φ(Ή ^~ . f, T)] _} = |A ^r 5 4- SA]# 4- (¾ ^" f . f) ^T fi0 - ti ^T S (H ^~ * _y T, ).

By using the inequality of equation (49) and taking into account the Lipshitz condition ll/( ^T i) -/i ^T )!i≤ ^λ !Ι ^τ ι™ ^T zlL the following relationship can be obtained:

V≤ t? (A ^T 5-4-S'A -r cS 4- sA/) ΰ , (51) and this will be negative, if and only if:

A ^T S ^? 4-SA 4-€λ ² Ι + s ^" -SS < 0. (52) which implies that the error dynamics are asymptotically stable. Taking into account the inequality in equation (52), and using the Schur complement, the inequality in equation (48) is obtained.

SIMULATION AND RESULTS

[0057] Simulations were performed using the model of the MD process to check the behavior of the response, and then to check the effect of the developed controller. The simulations were done with real process parameter values, and implemented over different locations among the process domain.

[0058] Model simulation. The model was simulated with the real parameter values listed in the table of FIG. 8. The temperature difference between the feed and the permeate subsystems 103 and 109 of FIG. 1 was analyzed, where the inputs are 60°C for the feed inlet and 20°C for the permeate inlet. The behavior of the D process with no imposed control will be unstable, and thus the temperature difference along the membrane boundaries will not be constant and does not fall within practical values. FIG, 7 is a plot of temperature with respect to membrane height illustrating the steady-state response of the open loop system. The temperature difference of the boundary layer temperature without controller function of the membrane module length is shown. The difference evolved with time for different locations along the membrane boundary layer, and thus the transmembrane heat flux was not fixed along the membrane length. In addition, if the temperature difference is not sufficient, no heat flux will be produced or transferred. The induced fluctuations in the temperature difference along the membrane boundaries can lead to instability in the generated heat flux and consume more energy. The instability of the produced heat flux can cause a fluctuation in the water production as well as significant problems for the process.

[0059] Controller simulation. The role of the nonlinear controller is to maintain the difference temperature along the membrane boundaries to be almost fixed around a desired reference. The value of this reference has been chosen to be within a range of 5°C to 15°C, following practical guidance for sufficient water production, in the case of the closed loop system, the difference of the boundary layer temperatures was kept fixed at the desired reference (e.g., equal to 5°C or 10°C in the simulations). At the same time, the values of the control inputs were physically reasonable and energy saving. Referring to FIG. 8, shown are surface plots illustrating an example of the behavior of the process at steady-state among the process domain, it shows the convergence of the output to the desired reference among ail the process domain.

[0080] Referring next to FIG. 9A, shown is an example of the temperature difference in different locations along the membrane module length. The first plot 903 shows the temperature difference with respect to a 5°C reference at the top of the membrane module, where the second plot 906 shows the temperature difference at the middle of the membrane module, and the third plot 909 shows the temperature difference at the bottom of the membrane module for the closed loop response. FIG. 9B shows the convergence of the temperature difference for the top (6%), middle (34%) and bottom (67%) locations on the membrane boundaries to the desired reference of 10°C. The convergence of the

temperature difference is accompanied with an acceptable evolution of the control inputs within the predefined range.

[0081] FIG. 0A shows the evolution of control law inputs with time for the case with the 5°C reference. The inputs are within economical energy supply ranges. The evolution of the error with time proves the asymptotic stability achieved in the Lyapunov function derivative. The evolution of the reference tracking error with time, where the response is asymptotically stable using the observer-based nonlinear Lyapunov controller, is shown in FIG. 10B. FIG. 10C depicts the time evolution of the control inputs for the case with the 10°C reference. Reaching the steady state phase, with no disturbance, the D process experiences a stable behavior among all its designated domains. FIG. 10D illustrates the effect of the controller on the process during the steady state phase with the 10°C reference.

[0062] Controller gain variation. The controller gain is a constant to smooth out the response of the controller, and reduce the oscillations that may occur in the transient response. The controller gain is a tradeoff between the convergence speed to steady state and the smoothness of the response. The effect of the controller gain on the control law is illustrated for the feed and permeate in FIGS. 1 1 A and 11 B, where it is noticeable that a low gain gives low oscillation control law. A low gain gives an oscillation free control law response. At the same time, reducing the gain forces the control law to be within the specified economical inputs range, without spikes in the response.

[0083] A Lyapunov based boundary controller for the nonlinear membrane distillation model has been presented. The model is based on a 2D advection diffusion equation and extended to be symmetric, and flexible in terms of the boundary conditions. Semi- discretization has been applied over the continuous PDE model, and a system of nonlinear ODEs has been obtained and written in a state-space representation, where it presents the model in a suitable form for control techniques. A good approximation for the continuous model can be granted, in which spatial step sizes have a significant role for stabilizing the numerical solution, and refining the model. The disclosed controller was applied over the state-space representation of the nonlinear model, where the output converges to the desired reference asymptotically. Certain limitations were put over the control law inputs to be within an economic energy saving ranges, and achieve asymptotic stability to the system. The effect of the controller gain was justified to be within small values to avoid oscillations in the responses. The Appendix provides additional description related to the control of the boundary layers temperature differences through the inlet temperatures and/or the inlet flow rates for the feed and permeate subsystems 103 and 109 of the D process 100 of FIG. 1. in some embodiments, the water production can be maximized (optimized) while decreasing the energy consumption of the D process 100. in other embodiments, the water production can be maintained at a desired level using the nonlinear Lyapunov controller.

EXPERIMENT AND RESULTS

[0064] A fully automated DCMD experimental lab scale set-up was locally designed and fabricated. FIG. 12 shows a schematic diagram of the experimental set-up. A flat sheet membrane module 1203 with an active surface area of 0.005 m ² , with flow channel dimensions of 0.1 m * 0.05 m * 0.002 m, was designed and fabricated using

polymethylmethacrylate (PMMA) material. A composite membrane with a

polytetrafluoroethylene (PTFE) active layer and a nonwoven polypropylene support layer was used for MD testing. Details of the membrane characteristics are presented in the table of FIG. 13A and the feed water (Red Sea water) characteristics are presented in the table of FIG. 13B. Red Sea water was preheated to desired temperatures and circulated through the feed side of the membrane module 1203, while de-ionized (Dl) water was circulated through the permeate side of the membrane module 1203 simultaneously in a countercurrent mode. Feed and permeate temperatures were controlled using thermoregu!ators (e.g., electric heater 1206 and electric chiller 1209). [0085] Fresh seawater was used as feed for each experiment. The permeate generated during the D process resulted in an increase in the volume of coolant and overflow 1212 was allowed to happen through the outlet of the permeate/coolant tank. The overflow permeate was collected in a separate container placed on a weighting balance. The increase in weight of the container was continuously monitored and recorded. Pressure (P), temperature (T) and flow rates (F) were continuously monitored at the inlet and outlet of the membrane module 1203 of both the feed and permeate sides using respective sensors connected to data acquisition unit equipped with a Lab View Software. Conductivities of both permeate and feed solutions were continuously monitored and measured using conductivity meters (e.g., Oakton Eutech instruments, Malaysia).

[0088] The influence of feed water temperature ranging from 40°C to 70°C on water vapor flux was studied by maintaining a constant temperature on the permeate side at 20 °C. Ail the data were recorded after reaching the steady state temperature. Similarly, the influence of feed water temperature ranging from 30 °C to 75°C on trans- membrane flux was also studied by ramping the feed solution temperature at a rate of 0.05°C per minute. During ramping experiments, the coolant temperature was also maintained at 20 °C. Feed and coolant flow rates were kept constant at 90 L/h and 60 L/h, respectively, for ail the experiments. The salt rejection was calculated by: SR = (1 - C _p /C _f ) * 100, where C _p and C _f are the salt concentrations of permeate and feed solutions, respectively. A salt rejection of 99.99% was observed in ail experiments. The water vapor flux was determined by J _v ~ mw/At, where mw is the weight of collected permeate at a particular time interval t, and A is the effective membrane area.

[0067] Model validation. The ADE model was simulated with a developed MATLAB® __ software and validated with experimental data through two sets of experiments. For validation purposes, the same MD parameter values and membrane specifications as well as environment conditions were applied. [0088] Steady-state validation. Referring to FIG, 14A, shown is the validation at the steady-state phase between the experimental data and the simulations for different feed inlet temperatures. The absolute error is also presented. The model acts best when compared to experiments at high feed inlet temperatures varying from 55°C to 7Q°C, while the absolute error is a bit larger for lower temperatures below 40 °C. This may be attributed to the fact that water vapor flux is higher at high feed temperatures, which reduces the errors in

measurements compared to the very low flux obtained at low feed inlet temperatures.

[0069] Transient validation. In the other set of experiments, the feed inlet temperature was increased in a ramp fashion from 30°C to 75 °C. The increment was 0.1 °C every 2 minutes, and the outlet feed temperature as well as the permeate flux were recorded. The comparison between the experiments and the model is shown in FIG. 14B. in this

observation, the transient response of the outlet feed temperature and the permeate flux showed high consistency with the experimental data. This contributes towards the validity of the dynamic model to be extended to the transient phase.

[0070] The MD process is highly sensitive to considerable variations in the membrane parameter values. For example, FIG. 14C shows the effect of the pore size variation on the membrane flux, where it is noticeable that the higher pore size leads to higher permeate flux. At the same time, the Liquid Entry Pressure (LEP) should not be exceeded to ensure no pore wetting is occurring inside the membrane pores. The LEP is equal to 15 PSi for the used membrane, while the actual pressure difference at the liquid-vapor interface AP is lower or equal to (7.64-8.33 PSi) for pore size equals to 0.5 μπ) and along temperature variation between (20°C ~80°C), respectively. The variation of the membrane parameter values is very likely to happen when manufacturing the membranes, and so it leads introduce some divergence in some cases.

[0071] Model analysis. The ADE model was applied inside the feed solution as well as the permeate solution. For simulations and up-scaling purposes a software in MATLAB© was developed to solve the dynamic model and obtain behavior plots for analysis and improvement aspects. In the simulations, real numeric values for the membrane parameters were used as listed in the tables in FIGS. 13A and 5. The simulations show the response in transient and steady-state phases. The feed and coolant temperatures were set at 60 °C and 20 °C, respectively.

[0072] FIG. 6A illustrates the evolution of the feed and permeate solution

temperatures at specific instant of time during the transient phase. The left plot of FIG. 16A shows the feed bulk temperature to be higher than the temperature at the membrane boundary layer from the feed side, and in a same manner the right plot of FIG. 6A shows that the temperature of the membrane boundary layer from the permeate has a higher temperature than permeate bulk, it also shows that the difference of temperature between the two sides of the membrane (the driving force) is higher at the top side of the module (feed inlet side) and decreases gradually towards the bottom side (feed outlet side), which is in agreement with the experimental data. FIG, 16B shows the evolution of temperature at the membrane interfaces, where the feed interface (curve 1403) ranges between 55 °C and 60 °C and the permeate interface (curve 1406) ranges between 20 °C and 25°C.

[0073] The flow rate of the feed and the permeate streams have a significant role in managing the heat transfer inside the D module1203 (FIG. 12). The flow rate controls the amount of the convective heat that transfers throughout the bulk solution, in which increasing its value makes the heat transfer via convection in the bulk solution more than the conduction through membrane pores. It is useful to have high convection effect, however continuous increasing in the flow rate is unprofitable, because the heat flux would saturate eventually, A smart compromise can take place between convection and conduction mechanisms. The relation between the feed flow rate and the temperature for different flow rates is illustrated in FIG. 16C. Temperature of water molecules near the inlet increases rapidly while increasing flow rate, and therefore, temperature polarization and residence time decrease, while energy consumed by the flow pumps increases. [0074] The temperature polarization was simulated to check the time evolution of its coefficient (TPC). The formula of the TPC is:

The value of TPC varies from zero for systems that suffer from large boundary layer resistances to unity for well-designed systems. FIG. 16D shows the evolution of TPC with time. During the transient phase of the process the amount of heat transfer through the membrane is low due to the starting of heat accumulation, and then it magnifies and increases after reaching steady-state phase.

[0075] A dynamic Lyapunov-based model was validated for control of the MD process. The mechanisms of heat transfer that occur in the feed and permeate solutions were modeled with a 2D Advection-Diffusion Equation (ADE) model. The model is able to describe the behavior of the MD process in a transient phase as well as in a steady-state phase. The mode! predictions were strongly correlated with experimental data in both steady state and transient phases, with model predictions being within ±5% of the experimentally obtained values.

EXTREIViUIVl SEEKING CONTROL (ESC)

[0076] The use of a perturbation-based extremum seeking control (PESC) scheme is proposed, with no loss of generality to all ESC schemes, to maximize the clean water production of a MD process with respect to the lowest energy consumption. The principal idea of the ESC is achieved by finding suitable operating set-points that are able to maximize or minimize a desired performance function. The PESC scheme can be used to simultaneously identify and regulate the MD process dynamic model to generate maximum permeate flux subjected to the minimum energy consumption of either the feed pump flow rate alone with constant permeate pump flow rate, or both the feed and the permeate pumps flow rates. The dynamic model for the D process based on the 2D ADE can be used. Moreover, the fact that the MD process is nonlinear with most states inaccessible, the nonlinear Lyapunov-based observer can be used to provide the PESC with an accurate estimate of the process states using the available measurements.

[0077] Perturbation based ESC. PESC is a scheme for the extremum seeking control (ESC). It has the advantage of performing continuous gradient type optimization with fast adaptation and easy implementation. PESC uses a dither signal (e.g., a sine wave) to perturb the performance function and estimates its gradient. Then, the gradient estimation is used to drive the process to a desired extremum using an integrator and the same dither signal. FIG. 7A illustrates an example of the components of a basic PESC scheme. A wash out filter is one of the components that is used to kill the DC component of the signal J.

Despite the fact that multiplying the signal y with a zero mean dither signal sin(a)t) performs the same role, the wash out filter is more effective and does not require the dither frequency ω to be relatively high. The low pass filter is not a necessary component for the PESC design, however it can be important for filtering out CQS(2&>£) signals.

[0078] The PESC scheme in FIG. 17B provides few design parameters (e.g., a, κ, ω, ω _ι and o _h ) to tune during the setting of the process. Some recommendations for tuning these parameters are include a smaller amplitude a can lead to smaller residual error at the minimum, but it may also lead to being stuck at a local minimum, and the converse can also be true. Moreover, the amplitude a jointly with the gain κ can control the convergence speed of the scheme. Wash out and low pass cut-off frequencies o½ and ω,, respectively, should be smaller than the dither signal frequency ω, so the filter kills the DC component of the signal J without distortion to the gradient estimation.

[0079] Problem formulation. Consider the MD process dynamic discrete model given by: where vecto T e R ^{2il' l} r is the siaie vector which represents the temperature distribution at each grid point) inside the feed and the permeate domain (T T _p ) for the time evolution as described with respect to equations (12)-(17), and the generated performance function:

for optimization. The performance function, ), is a function of the temperature distribution inside the designated process domain T via the permeate vapor flux J(T), especially the ones associated with boundary layers of the membrane, and u e [u _f , u _p ] ^T , which represents the pumps flow rates for either the feed container or both the feed and the permeate containers, p is the water density, and it is used to fix the unit consistency for the performance function, μ is a regularization coefficient that is used to adjust the weight of each term in the performance function. During the process simulations, it is possible to shed more light on certain cases where one term is more dominant over the other. iJ(T) j is the absolute value of the permeate vapor flux as in the Knudsen diffusion model for the mass transfer coefficient given above, where the negative sign indicates the maximization.

[0080] Single parameter PESO . Consider optimizing the performance function in equation (55) to be subject only to the feed pump flow rate. Consequently, permeate pump flow rate u _p is kept constant, while u = u _f is the control. Note that in practice, the feed pump flow rate has significant influence on the process performance more than the permeate pump flow rate.

[0081] The PESC scheme is applied to equation (55), where the permeate flux function is evaluated simultaneously with the evaluation of the temperature distribution (T) using the 2D ADE model. Afterwards, the DC content inside the performance function in (55) is eliminated using the wash out high pass filter, as shown in:

^■ ί? ) ( ÷ ¾) ~

(56) Then it is multiplied with the dither signal a sin(cot) that has a small amplitude a, and an integrator with gain κ to extract the gradient direction. u ^~ a Ksin ( cot ) u ^{■■■■■} a sin ( mi ) . (57)

The adaptation signal, u, shifts the dither towards the gradient direction, as described in:

It is recommended that the cut-off frequency of the wash out high pass filter ω _Λ , to be smaller than the frequency of the dither signal ω, and the amplitude of the dither signal a and the gain κ be small as well. FIG. 17B illustrates the principal steps that the PESC follows to update the control.

[0082] Multiple-parameter PESC. To further broaden the optimization analysis on the D process, consider optimizing the performance function in (55) to be subject to both the feed and the permeate pumps flow rates to be control inputs. This would considerably highlight on the impact of the feed pump flow rate and the permeate pump flow rate on the process production, and the interference between them.

[0083] The optimization procedure is similar to the single parameter case, however the multiple-parameter PESC considers each control input in a separate loop with a specific dither signal that has some phase shifts (e.g. 0 and 772). These phase shifts are introduced to compensate for the delays provided by the plant dynamics and the performance function J(T, u) in equations (54) and (55), respectively. The design parameters are set for each input independently, and then the whole process is tuned out to have an optimized behavior. The following equations show the dynamic representation for the feed and the permeate pumps flow rates. f = a-2 i¾sw(<%f— <j>2 ) (¾ J) (59) Up ~~ φ i ) ( η i J) (60)

(61) i|i ^ ^~~ mhi + ⁽ *>hiJ (62) Afterwards, each control inputs will consist of the adaptation signal in addition to the dither signal.

U = ^f + sinim t} (63) Up = itp + aisi i&it) (64)

FIG. 1 7C illustrates a diagram of the main components for the multiple-parameter PESC.

[0084] Anti-windup PESC. The presence of an integrator in any PESC causes, in some occasions, the control inputs to wander outside their acceptable and feasible domain. This may facilitate the actual inputs to saturate on the actuators (e.g., pumps) domain boundary. This phenomenon is similar to the integrator wind-up in classic control. Consider an anti- windup technique to overcome input saturation phenomenon. This technique is preferable over dealing with the saturation as constant, since the process is dynamic and more likely to have input saturation. This technique allows the imposition of certain constraints on the control inputs to bound them with physical limits. Technically, it drives the actuator back to the unsaturated mode when input saturation is about to happen by using a penalty-like scheme to prevent or overcome saturation.

[0085] Combine the anti-windup technique along with the PESC scheme, and define a saturation function as in:

f &max tf

sat (tt, u _s ) ~ - I u if ttmin ^li U ax , (65)

1 H-min tf It °E Itfjijfj

where u _s = [u _max , u _min ] ^r and u _max and u _mtn are the upper input limit and the lower input limit, respectively. K _aw is a design parameter that adds more significance to the saturation function weight. The saturation function is placed before the process model, and defines an appropriate input v for a feedback. FIG. 17D shows a diagram for a single parameter PESC with anti- windup technique. [0086] Similar to the previous cases, the adaptation signal is evaluated, however in this case, it includes the effect of the anti-windup. Following equations present the adaptation signal.

= ai qsin(€M)u— α$κ$ίη(ωή

+ Kaw (sat (i _s ) ~ //), (66)

The same procedures are applicable to multiple-parameter PESC, however the anti-windup model is added to each control input channel to prevent input saturation.

OBSERVER BASED PESC

[0087] The MD process model is observable, and the present nonlinearity is Lipshitz continuous. For these reasons, an observer can be used to provide an estimate of the MD process states and feed them to the PESC. The observer is nonlinear due to the presence of nonlinearity in the process model of equation (54).

[0088] Nonlinear Lyapunov-based observer design. The advantage of having a good model for the present nonlinearity enables consideration of a nonlinear observer that can be modified with an algorithm to find an adequate observer gain. Therefore, the D process nonlinear model in equation (54) can be extended and adapted, where the stability proof of the observer was done for a class of nonlinear systems given by:

f _^ AT + Bu+fiT). (68) where A is the dynamics matrix, B is the input matrix, f(T) is the nonlinearity vector in the model and T is the state vector. For the observer design, consider the nonlinear Lyapunov- based observer in the following form:

T(x> ) = A(«)T(A ^' _? ) -r /(T(*,z,/)) + L0 ^{» ~} CT) . (69) where T e Μ ^2η2χ1 is the estimated state vector, and L is an observer gain matrix. Then, the observer error dynamic equation is given by:

T = (A ····· LC)f +/(T) -/(f) , (70) where T = T - T is the state estimation error. Considering Proposition 2, the foiiowing provides sufficient conditions for exponential stability of the state observer error of equation (70).

[0089] It can be shown that if X and Y are real vectors of the same dimension, then the following inequality holds:

— 4

Define A ₀ = A( ) - LC, which is valid since A( ) is constant for the observer use, and a Lyapunov function as:

V - f ⁷ ΡΊ (72) for the stability, and its derivative as:

V - ~T ^r [Aj - PA ₀ ]T - 2T ^r P[/(T) - /(f)], (73) where P e M ^{2n' x2n^} is any arbitrary symmetric positive definite matrix. Then, the following is valid since (T) is Lipschitz.

2 ^r Pi/ ^' (T) - /(f )] < 2|PT|y|f I . ₍₇ 4)

Using the inequality of equation (71), define:

Then,

2y|FT[ |f| < r T ^T P T ⁷ + T ^r f, (76) V < Ψ [A| + PAo + Γ P - /]T.. (77) Note that if [A _Q I + Pi ₀ + γ ² ΡΡ + /] is negative definite, then there exists a P ₂ such that:

[Aj I 4- i¾A ₀ -f fPi + 1} < (78)

Let

Ao P + PAo -f†PP - = --£/. (79) Then:

V < -εΎ ^τ Ύ < 0 so T - 0 exponentially. (80) (A0) - LC) P -f P(A («) - LC + PP + /(ε -f 1 ) = 0 (81)

By setting the terms as described in Proposition 2, a standard algebraic Ricatti equation is obtained that has a solution if P is a symmetric and positive definite matrix.

[0090] FIG. 8A shows the fast convergence of the estimated feed (T _f ) and permeate (7 _p ) states using the nonlinear Lyapunov-based observer for different locations in the MD discretization grid. FIG. 18A also shows how the nonlinear observer is robust to any sudden changes in the process. The fast convergence of the observer is guaranteed by taking smaller sampling time, in the same context, the error analysis depicted in FIG. 18B illustrates the exponential convergence of the error (AT) to zero.

[0091] Observer based PESC structure. Following the design of the nonlinear

Lyapunov-based observer, the temperature distribution ail along the designated process domain can be available at each instant of time. Proceeding further, these estimates can be provided to the PESC to optimize the performance function in equation (55), without the need to directly access the inaccessible process states. This is valid, when the observer convergence is fast.

[0092] Considering the single parameter PESC, the structure of the performance function would be given by:

and the PESC structure would depend only on the estimate of the temperature distribution (T). FIG. 19 illustrates an example of the structure of the observer based single parameter PESC.

SIMULATION AND RESULTS

[0093] Analysis for the observer based PESC scheme will now be discussed through numerical simulations using the MD process dynamic model in equation (54) for each input case. The PESC scheme was equipped with anti-windup technique to prevent input saturation. All fabrication parameters of the employed membrane were realistic, and are listed in the table of FIG. 20. For the single parameter case, we consider the feed pump flow- rate (u = U ) as the only control input, while keeping the permeate pump flow rate constant at u _p = 0.2 (m/s). From a practical perspective, u _f should be bounded with minimum and maximum values to reflect the pump's physical limitations and energy management, as given by 0.1 (m/s) < u _f < 0.4 (m/s),

[0094] The simulation results of the single parameter PESC with anti-windup technique is depicted in FIG. 21A, where it illustrates the time evolution of the performance function, the permeate flux and the optimized feed pump flow rate (u _f ). It can be seen that the permeate flux reaches its maximum value (about 12 kg/m ² hl, and the control input oscillates in a small neighborhood around the optimal value, which is less than its maximum bound.

[0095] Referring next to FIG. 21 B, shown are the simulation results when a sudden disturbance is introduced to the D process to assess the sensitivity as well as the adaptability of the PESC scheme. The disturbance happens at the middle of the simulation by lowering the feed inlet temperature from 60 °C to 30°C at t = 50 s. FIG. 21 B shows the response of the control input to the disturbance. Note that the generated permeate flux is more significant at high temperature, and this may explain the low generation of permeate flux at 30 °C.

[0096] For the multiple-parameter case, the performance function was optimized considering the feed and the permeate pumps flow rates as control inputs, where u =

[u _r , . FIG. 21 C demonstrates the time evolution of the performance function, the permeate flux and the optimized feed and permeate control inputs. FIG. 21 C also shows the rapid convergence of the permeate flux to its maximum value, and the saturation of the control inputs. It can be seen from the simulations results, that the permeate flux in the multiple-parameter case reaches a similar level as the single parameter case (around about 12 kg/m ² h), but with higher pump flow rates. This may be attributed to the additional degrees of freedom thai are provided to the process, by assigning a second control input, where the newly added design parameters can be elaborately tuned, and therefore, contribute toward fast convergence for the performance function.

[0097] An observer based PESC for a D process has been examined. The process can be modeled with a dynamic 2D ADE model. It is possible for the controller to maximize the permeate vapor flux while at the same time to minimize the consumed energy via the feed and permeate pumps flow rates. For this, a nonlinear Lyapunov-based observer was used to provide an estimate for the temperature distribution among the domain of the MD process and feed the estimates to the PESC. Two sets of control inputs were considered for the PESC scheme: single control input of the feed pump flow rate, and multiple control inputs of both feed and permeate pumps flow rates. They were provided with an anti-windup technique to prevent actual input saturation. Simulations show the time evolution of the control inputs and support the feasibility of their energy consumption.

MULTIVARIABLE ESC OF SOLAR POWERED MD PROCESS

[0098] A real time multivariab!e optimal controller is disclosed for a solar powered direct contact membrane distillation (DCMD) process using Newton-based extremum seeking control (ESC). This controller takes into account the varying feed inlet temperature, due to the dynamic nature of the sun. Moreover, since ESC is non-model based, the controller design is robust against plant-model mismatch. For this work, a muitivariable Newton-based ESC is described for the optimization of a solar powered MD water desalination plant. The dynamic model of the heat and mass transfer in the D process along with the solar powered MD optimal control problem formulation will be presented.

[0099] First consider the MD module illustrated in FIG. 22. in this configuration, hot water passes along the hydrophobic membrane on one side, which is called the feed side, and cold water flows in the counter direction along the other side, which is called the permeate side. Water vapor is driven from the feed side across the membrane and into the permeate side by the induced vapor pressure difference. Both heat and mass transfer processes occur simultaneously as water evaporates at the feed-membrane interface and condenses at the permeate-membrane interface. As a result, the temperature at the membrane boundary layers differs from the bulk temperature of the feed and permeate streams. The D module can be divided into N control volume cells, where individual cells can have a uniform bulk temperatures T _bfn , T _bprt ) except at the boundary layers where the temperatures are (T _mfri , T _mp , as shown in FIG. 22 for the n ^fh ceil.

[0100] Mass transfer in DCMD. The transport phenomena is described by the classic gas permeation and heat transfer theories. The mass flux (J) in DC D is related to the saturated vapor pressure difference across the membrane (AP) through the membrane mass transfer coefficient (B _m ) as follows

The saturated vapor pressure of pure water (Pw ^at [T]) as a function of temperature is given by the Antoine equation:

Dissolved salt in the feed stream reduces the saturated vapor pressure, therefore to compensate for this the following relation was proposed in:

P _mf ^™ (1 - arj a0 )(l— Q& tfaCl ~ l0x _{N aCl} ) P ^* ? ¹ [1 inf _. (85) where x _Naa is the mole fraction of NaCi in the feed stream. However, permeate and the saturated vapor at the membrane-permeate interface is pure water. Thus, P _mp = .

[0101] Heat Transfer in DCMD. Heat transfer in DCMD is a spatially and temporally distributed phenomena, A dynamical model, based on an electrical analogy to thermal systems, can be used for both transient and steady-state response. The resulting model can be a system of differential-algebraic equations (DAE) or a descriptor system. The schematic diagram for the electrical analog circuit for the DCMD module is shown in FIG. 23. By applying Kirchoff s electrical laws (KCL and KVL) to the analog circuit for the DCMD, the following equations can be obtained. On the feed side, the rate of change of the heat transfer rate from the n-1 ceil to the n ^th cell is proportional to the temperature difference between them. Taking into consideration the series impedance Z„, this can be expressed as:

^™ ¾ ^™ ^ Jn ¹ bfn-i ~ j % - JT bt _u . (86)

Using Kirchoff's current iaw at the n ^th feed node, it foiiows that the rate of change for the bulk feed temperature (¾ _η ) is:

^ ^" ά (i ^{+ JnAmCp} ) ^Th{*

Equation (88) describes the dynamics of the heat transfer rates into and out of the n ^tn feed cell ((?;·„ and , respectively).

[0102] On the permeate side, the rate of change of the heat transfer rate (Q _Pn ) is:

And the dynamics of the bulk permeate temperature (T _b? is?

at - ~ bp p - bp

The coupling between the feed and the permeate dynamics in the rf ceil can be written in residue form as:

The outlet temperatures at the terminal ceils of the feed and permeate analog are also given by the algebraic equations, where are respectively:

0— 2jf _ou¾ — 2 _Pin - - ¾ _er∞ Qf _a .i? (91)

0 = ip _out — If _in + ·¾· _;¾ί;η Q : · ₉₂ )

[0103] The heat and mass transfer equations (88)-(92) represent a nonlinear differential-algebraic system. When considering N number of interconnected ceils, the resultant equations can be expressed as a nonlinear descriptor system of the form:

EX it) ^ F(X(t), u t)) (93)

J ^ gi ) (94) where X E R ^6JV+4 represents the differential and algebraic states X refers to the time derivative of the state vector, E ε jt ^{6W +4x6iV+4} is a singular matrix, F(X(t), u(t)) is a nonlinear function of the states and input, (t) = [Mf _ipl M _Pin ] represents the inlet feed and permeate mass flow rates respectively, the distilled water flux J is also given by a nonlinear function of the state g(X). Wth the dynamical model of DCMD described, next consider the complete setup of the solar powered DCMD water desalination system.

[0104] Optima! control problem of solar powered DCMD. Referring to FIG. 24, shown is an example of a solar powered MD process. The solar powered DCMD water desalination setup comprises a solar thermal collector 2403 connected a storage tank 2406 and a heat exchanger 2409. The feed water is circulated from a feed tank 2412 through the heat exchanger 2409 to the DCMD module 2415 and back to the feed tank 2412. The permeate is pumped from a fresh water tank 2418 into the DCMD module 2415, where it collects the distilled water flux and flows back into the fresh water tank 2418. The control inputs are the mass flow rates of feed and permeate sides. The feed inlet temperature (7 _f . _n ) will vary according to the solar radiation throughout the day, which affects the distilled water flux. An optimal control strategy can be used to ensure the efficiency of the process under disturbed operating conditions. One way to optimize the process is by maximizing the following objective function Y, which minimizes the feed/permeate inlet mass flow rates and maximizes the distilled water flux (/): max J - max ( χ,7 - (ο÷2-Μ¾ _η 4- . (95)

[0105] FIG. 25A shows the response of the DCMD process under various feed inlet temperatures and feed inlet mass flow rates. In this simulation, the permeate inlet mass flow rate was kept constant at 256 Kg/hr, while the feed inlet mass flow rate was increased from 90 Kg/hr to 375 Kg/hr. This was repeated for 5 feed inlet temperatures. As demonstrated in FIG. 25A, increasing the feed inlet temperature increased the distilled water flux. While increasing the feed inlet mass flow rate increases the flux, the flux reached saturation values at high feed inlet mass flow rates. As shown in FIG. 25B, the objective function of equation (95) has a maximum that varies for different feed inlet temperature. Similar characteristics are exhibited for the objective function with respect to the inlet permeate mass flow rate. Therefore, the controller can automatically track the peak of the objective function by manipulating the inlet feed and permeate mass flow rates,

iViULTi VARIABLE ESC FOR REAL TIME OPTIM IZATION

[0106] Extremum seeking is a real time optimization method, which is particularly useful to keep the reference-to-output map at its extremum (maximum or minimum). This extremum might shift with the process parameters or under disturbance influence and it's the job of the ESC to track the extremum as if shifts. Newton-based multivariable ESC can rely on sinusoidal perturbation signals to estimate the gradient and the Hessian of the objective function by using low and high pass filters. It has been shown that the objective function has a (7 . _n ) dependent peak, therefore it can be assumed that there exists u ^* = such that: ^■ A ) (96)

chiwhere H = H ^'r .

[0107] Suppose that the control input is given by u = β(Χ, Θ), which is a smooth function of the states and the parameters θ = [θ _1ι θ ₂ ] ^Γ . Then, the dosed loop system is given by EX = F(X, (X, 0)) . The objective of the ESC is to maximize the value of Y. The block diagram of the Newton-based ESC is shown in FIG. 28, where K _n is a positive diagonal matrix, and the perturbation matrices (L(t) e M ^2xl , Ω(£) e M ^2x2 , and (t) e M ^2xi ) are defined as:

0

Lit ) — sin^ _! t)<— sm{u <yt } (98)

4

Qu ( s (a¾t) sin(u¾t)) , i φ k (101) i * (102) rn (t)— [<¾ sm(wi*) _i Q2 sin(w2i) ^" P , where ω _έ ≠ s _k and such that ωι/ω is a rational number, and c _1; c ₂ are real positive numbers. The cut-off frequencies for the low and high pass filters, ω ₍ and ω _η respectively, can be designed appropriately according to:

(103) ( 04)

This design derives an estimate of the gradient vector G and the Hessian by adding the perturbation signal φ(ΐ) to the estimated optimal input Θ. The estimate of the Hessian matrix is inverted by the dynamical system given by:

where Γ = H ^"1 , to avoid difficulties of algebraically inverting H when it is dose to singular. To summarize the dynamical system shown in FIG. 26, the following equations are presented:

P / . / . P DT where

SIMULATION AND RESULTS

[0108] Simulations were carried out in order to demonstrate the effectiveness of the ES controller to optimize the solar powered DCMD water desalination process. For this purpose, realistic membrane parameters were used as listed in the table of FIG. 27. The design parameters for the Newton-based ESC were: c = c ₂ = 0.15, ω ₁ = 3 rad/s, ω ₂ = 2 rad/s, ω, = 0.1 rad/s, ω _κ = 4 rad/s, K _n = diag([0.05, 0.05]), «) _r = 0.0008.

[0109] The feed inlet temperature profile depicted in FIG. 28A was designed to reflect a realistic behavior. The permeate inlet temperature was kept constant throughout the simulation at 20°C. As the feed inlet temperature increases, the distilled water flux (/) is optimized taking into account the feed and permeate inlet mass flow rates, as shown in FIG. 28A. FIG. 28B shows the estimated optimal inlet mass flow rages for the feed and permeate sides along the objective function value. The ESC dynamics takes about 50 seconds to converge to the optimal estimated input Θ as shown in FIG. 28B. After that, the controller tracks the peak of the objective function with smooth transitions as it varies with the feed inlet temperature. As the temperature drops, the ESC quickly adjusts the inputs.

[0110] The optimization of the solar powered DC D water desalination process was presented. The optimization problem addressed maximizing the distilled water flux while minimizing the feed and permeate inlet mass flow rates, which are highly related to the energy being used to pump water. The multivariabie Newton-based ES controller was used for real time optimization. A dynamic model of the process was used to carry out the simulation. The simulation results show that the ESC was able to track the peak of the objective function as if varied with the feed inlet temperature.

[0111] With reference now to FIG. 29, shown is a schematic block diagram of an example of processing circuitry 1100 that may be used to implement various portions of the control of the membrane distillation (MD) process 100 of FIG. 1 in accordance with various embodiments of the present disclosure. The processing circuitry 1100 includes at least one processor circuit, for example, having a processor 1103 and a memory 1 06, both of which are coupled to a local interface 109. To this end, the processing circuitry 100 may be implemented using one or more circuits, one or more microprocessors, microcontrollers, application specific integrated circuits, dedicated hardware, digital signal processors, microcomputers, central processing units, field programmable gate arrays, programmable logic devices, state machines, or any combination thereof. The local interface 1 09 may comprise, for example, a data bus with an accompanying address/control bus or other bus structure as can be appreciated. The processing circuitry 1100 can include a display for rendering of generated graphics such as, e.g., a user interface and an input interface such, e.g., a keypad or touch screen to allow for user input. In addition, the processing circuitry 100 can include communication interfaces (not shown) that allow the processing circuitry 1100 to communicatively couple with other communication devices. The communication interfaces may include one or more wireless connecfion(s) such as, e.g., Bluetooth or other radio frequency (RF) connection and/or one or more wired connection(s). [0112] Stored in the memory 1 106 are both data and several components that are executable by the processor 1 103. In particular, stored in the memory 108 and executable by the processor 1 103 are MD process application(s) 1 115, an operating system 1 118, and/or other applications 1121. MD process applications 1 115 can include applications that support, e.g., controllers for control of the operation of the MD process 100 and/or observers for estimation of states of the MD process 100. It is understood that there may be other applications that are stored in the memory 1106 and are executable by the processor 103 as can be appreciated. VVhere any component discussed herein is implemented in the form of software, any one of a number of programming languages may be employed such as, for example, C, C++, C#, Objective C, Java ^® , JavaScript ^® , Perl, PHP, Visual Basic ^® , Python ^® Ruby, Delphi ^® , Flash ^® , LabVIEW ^® or other programming languages.

[0113] A number of software components are stored in the memory 106 and are executable by the processor 1 103. in this respect, the term "executable" means a program file that is in a form that can ultimately be run by the processor 1 103. Examples of executable programs may be, for example, a compiled program that can be translated into machine code in a format that can be loaded into a random access portion of the memory 106 and run by the processor 1103, source code that may be expressed in proper format such as object code that is capable of being loaded into a random access portion of the memory 1106 and executed by the processor 1 03, or source code that may be interpreted by another executable program to generate instructions in a random access portion of the memory 1106 to be executed by the processor 1103, etc. An executable program may be stored in any portion or component of the memory 106 including, for example, random access memory (RAM), read-only memory (ROM), hard drive, solid-state drive, USB flash drive, memory card, optical disc such as compact disc (CD) or digital versatile disc (DVD), floppy disk, magnetic tape, or other memory components,

[0114] The memory 1 106 is defined herein as including both volatile and nonvolatile memory and data storage components. Volatile components are those that do not retain data values upon loss of power. Nonvolatile components are those that retain data upon a loss of power. Thus, the memory 1 108 may comprise, for example, random access memory (RAM), read-only memory (ROM), hard disk drives, solid-state drives, USB flash drives, memory cards accessed via a memory card reader, floppy disks accessed via an associated floppy disk drive, optical discs accessed via an optical disc drive, magnetic tapes accessed via an appropriate tape drive, and/or other memory components, or a combination of any two or more of these memory components, in addition, the RAM may comprise, for example, static random access memory (SRAM), dynamic random access memory (DRAM), or magnetic random access memory ( RAM) and other such devices. The ROM may comprise, for example, a programmable read-only memory (PROM), an erasable

programmable read-only memory (EPROM), an electrically erasable programmable readonly memory (EEPROM), or other like memory device.

[0115] Also, the processor 1 103 may represent multiple processors 1 103 and the memory 1108 may represent multiple memories 1106 that operate in parallel processing circuits, respectively, in such a case, the local interface 109 may be an appropriate network that facilitates communication between any two of the multiple processors 1103, between any processor 1103 and any of the memories 1108, or between any two of the memories 1 106, etc. The local interface 1109 may comprise additional systems designed to coordinate this communication, including, for example, performing load balancing. The processor 1103 may be of electrical or of some other available construction.

[0116] Although the MD process application(s) 11 15, the operating system 1 118, application(s) 1121 , and other various systems described herein may be embodied in software or code executed by general purpose hardware as discussed above, as an alternative the same may also be embodied in dedicated hardware or a combination of software/general purpose hardware and dedicated hardware. If embodied in dedicated hardware, each can be implemented as a circuit or state machine that employs any one of or a combination of a number of technologies. These technologies may include, but are not limited to, discrete logic circuits having logic gates for implementing various logic functions upon an application of one or more data signals, application specific integrated circuits having appropriate logic gates, or other components, etc. Such technologies are generally well known by those skilled in the art and, consequently, are not described in detail herein,

[0117] Also, any logic or application described herein, including the MD process application(s) 11 15 and/or application(s) 1121 , that comprises software or code can be embodied in any non-transitory computer-readable medium for use by or in connection with an instruction execution system such as, for example, a processor 1 103 in a computer system or other system. In this sense, the logic may comprise, for example, statements including instructions and declarations that can be fetched from the computer-readable medium and executed by the instruction execution system. In the context of the present disclosure, a "computer-readable medium" can be any medium that can contain, store, or maintain the logic or application described herein for use by or in connection with the instruction execution system. The computer-readable medium can comprise any one of many physical media such as, for example, magnetic, optical, or semiconductor media. More specific examples of a suitable computer-readable medium would include, but are not limited to, magnetic tapes, magnetic floppy diskettes, magnetic hard drives, memory cards, solid-state drives, USB flash drives, or optical discs. Also, the computer-readable medium may be a random access memory (RAM) including, for example, static random access memory (SRAM) and dynamic random access memory (DRAM), or magnetic random access memory (MRA ). In addition, the computer-readable medium may be a read-only memory (ROM), a programmable read-only memory (PROM), an erasable programmable read-only memory (EPROM), an electrically erasable programmable read-only memory (EEPROM), or other type of memory device.

[0118] It should be emphasized that the above-described embodiments of the present disclosure are merely possible examples of implementations set forth for a clear understanding of the principles of the disclosure. Many variations and modifications may be made to the above-described embodiment(s) without departing substantially from the spirit and principles of the disclosure. All such modifications and variations are intended to be included herein within the scope of this disclosure and protected by the following claims.

[0119] It should be noted that ratios, concentrations, amounts, and other numerical data may be expressed herein in a range format. It is to be understood that such a range format is used for convenience and brevity, and thus, should be interpreted in a flexible manner to include not only the numerical values explicitly recited as the limits of the range, but also to include all the individual numerical values or sub-ranges encompassed within that range as if each numerical value and sub-range is explicitly recited. To illustrate, a concentration range of "about 0.1 % to about 5%" should be interpreted to include not only the explicitly recited concentration of about 0.1 wt% to about 5 wt%, but also include individual concentrations (e.g., 1 %, 2%, 3%, and 4%) and the sub-ranges (e.g., 0.5%, 1.1 %, 2.2%, 3.3%, and 4.4%) within the indicated range. The term "about" can include traditional rounding according to significant figures of numerical values. In addition, the phrase "about 'x' to y includes "about 'x' to about 'y"\