Exploiting meteorological forecasts for the optimal operation of algal ponds

Journal of Process Control 55 (2017) 55–65

Contents lists available at ScienceDirect

Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont

Exploiting meteorological forecasts for the optimal operation of algal ponds Riccardo De-Luca a , Fabrizio Bezzo a , Quentin Béchet b , Olivier Bernard b,∗ a CAPE-Lab Computer-Aided Process Engineering Laboratory, Department of Industrial Engineering, University of Padova, via Marzolo 9, I-35131 Padova, Italy b Université Côte d’Azur, Inria, Biocore, BP93, 06902 Sophia-Antipolis Cedex, France

a r t i c l e

i n f o

Article history: Received 30 October 2012 Received in revised form 22 August 2013 Accepted 23 August 2013 Available online 9 May 2017 Keywords: Optimization problems Optimal control Model Predictive Control Microalgae Temperature Growth model Biofuels

a b s t r a c t Biofuel production from microalgae requires optimizing the operation of cultivation systems (i.e. outdoor raceway ponds) for this process to be economically sustainable. Controlling algal ponds is complex as the cultivation systems are exposed to fluctuating conditions. The strategy investigated in this study uses weather forecasts coupled to a predictive model of algal productivity to optimize pond operation. Productivity was optimized by dynamically controlling rates of fresh medium injection and culture removal into and from the pond. This optimization strategy when applied to a cultivation plant in Nice, South of France, increases the productivity by 2.13 compared to the reference case where the pond depth and dilution rate were kept constant over time. The underlying Model Predictive Control consists of playing with raceway pond thermal inertia and supply of fresh water to reach rapidly optimal temperature, and then keep a balance between photosynthesis and respiration in the darkest layers of the raceway pond. The meteorological inaccuracy for forecasts beyond 24 h was compensated by frequent updates of the optimal control problem. Finally, this scheme turned out to be robust to inaccurate weather forecasts, and the net productivity value reached was close to the productivity obtained for perfectly known meteorology. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Microalgae have emerged during the last decade as one of the most promising new technologies for providing innovative molecules for the cosmetic and pharmaceutical industries, and as a source of proteins for animal and human nutrition [1,2]. At larger time horizon, microalgae are expected to contribute to fossil carbon replacement with renewable carbon, especially for supplying green chemistry and liquid biofuel in the transport sector [3]. The great interest in this technology is not only the substantial higher productivity compared to terrestrial plants, but also the possibility of coupling the microalgal production process to industrial CO2 mitigation and wastewater treatment to finally recycle carbon, nitrogen and phosphorus. Nevertheless, further developments and optimization steps are required to reduce both the cost and the environmental footprint of these processes. One way to increase the profitability of the system is to develop control strategies to maximize productivity

∗ Corresponding author. E-mail address: [email protected] (O. Bernard). http://dx.doi.org/10.1016/j.jprocont.2017.03.010 0959-1524/© 2017 Elsevier Ltd. All rights reserved.

while better managing the nutrient and water use. Controlling an outdoor algal production system is complex because the state variables of the cultivation system (e.g. biomass concentration, pond temperature, etc.) are continuously driven by the meteorological conditions (solar irradiance, air temperature, etc.). A first approach was derived by [4] who used the Pontryagin’s maximum to identify the set of necessary conditions to be satisfied in order to guarantee optimal control in a periodically varying environment. Nevertheless, the simplifications required by this approach, which assumed a known and periodic meteorology, reduced the applicability of the proposed optimal strategy in the case of algal outdoor cultivation. Model Predictive Control (MPC) [5] has proven to be very efficient to manage situations with complex modelling, where classical model based control cannot be easily derived. MPC seems especially relevant when meteorology plays a key role, such as energy efficient building climate control [6] or management of distributed power system with wind turbines [7]. MPC was already applied to microalgae [8–10], but never to manage weather forecasts in the control strategy. Alternatively, we propose a MPC strategy accounting for the knowledge of future weather conditions. This technique is based on mathematical models able to predict algal productivity from

56

R. De-Luca et al. / Journal of Process Control 55 (2017) 55–65

the knowledge of weather conditions [11,12]. In practice, this optimization task consists of determining the optimal rates at which fresh medium has to be injected into or extracted from a mixed algal pond (also named “raceway”) based on weather forecasts. The first objective of this study was to propose an improvement to the standard management approach, and assess the productivity gain of this optimized approach. The secondary objective was to understand and analyze the rational behind the predictive controller. We first detail the model consisting in a coupled biologyphysical model. Then we define the model predictive controller and we detail the numerical approach developed for solving the problem. We first consider the case where weather forecasts are perfect, and we analyse the logic behind the control action. Finally, we consider the more realistic case where meteorology becomes uncertain after 24 h.

2. The microalgal raceway pond productivity model The model couples two submodels: a physical submodel predicting the temperature in the raceway pond and a biological submodel computing microalgal production. Both submodels are driven by meteorological conditions determining the fluxes of heat and photons through the system. The model has a triangular structure since the biological submodel does not influence the thermal dynamics of the physical submodel. The model is a simplified version of the model in [13]. However in [13], the pond depth is maintained at a constant level, while here we exploit the idea of changing the depth to modify the thermal inertia of the system. The scheme of the model structure (coupling biological and thermal equations) and the objective function is shown in Fig. 1. The model equations are detailed in the following paragraphs.

2.1. Deriving the biological model from mass balance The biological model describes the dynamics of a microalgal biomass (which concentration is x(t) in kg m−3 ). Microalgae are photosynthetic microorganisms which capture CO2 in the presence of light and incorporate this carbon into their biomass, while they continuously loose a fraction of their carbon by respiration. Both activities are temperature dependent. The open algal raceway pond of volume V (m3 ) and depth lp (t) (m) is supposed to be perfectly mixed [14,15]. The total biomass x(t)V(t) in the open raceway pond varies over time according to the following equation: d(x(t)V (t)) = −x(t)qout (t) + G( · )V (t) − R( · )V (t), dt

(1)

where t is the time variable (s), G(·) is the depth-averaged specific growth rate (kg m−3 s−1 ) and R(·) is the specific respiration rate (respiration causes carbon loss) (kg m−3 s−1 ). The fresh medium is injected at the rate qin (t) (m3 s−1 ) whereas the culture is extracted from the raceway pond at the rate qout (t) (m3 s−1 ). The raceway pond volume varies over time according to the following equation: qin (t) − qout (t) + vr (t)S − me (t)S dV (t) = , w dt

(2)

where S is the raceway pond surface area (m2 ), w is the water density (kg m−3 ), vr (t) is the rainwater flow (m s−1 ), and me (t) is the evaporation mass flux (kg m−2 s−1 ). The average specific growth rate G(·) in Eq.(1) depends on the biomass concentration x(t), the raceway pond temperature Tp (t), and the solar irradiance Hs (t) (W m−2 ). The resulting function G(t, x(t), Hs (t), Tp (t)) was expressed as [16,17]: G(t, x(t), Hs (t), Tp (t)) =

1 lp (t)



lp (t)

m x(t) 0

H Hs (t)e−x(t)z KI + H Hs (t)e−x(t)z

dz (3)

where m is the maximum specific growth rate (s−1 ),  is the extinction coefficient (set at 120 m2 kg−1 ), due to light absorption and scattering, H is the photosynthetically active radiation (PAR) fraction of solar light (set at 0.47 [17]), z is the local depth (m) and KI is the half-saturation parameter (W kg−1 ). The impact of photoinhibition on microalgal growth was not explicitly included in this study. Indeed, as suggested by [11], a Monod kinetics can efficiently represent algal growth at high biomass density. This is explained by the fact that only a small fraction of cells are photo-inhibited in the dense cultures, leading to an average behaviour of Monod type. The specific respiration rate R(·) in Eq.(1) depends on raceway pond temperature Tp (t) and biomass concentration x(t) through the following equation: R(t, x(t), Tp (t)) = r (Tp )x(t),

(4)

where r is the respiration coefficient (s−1 ). The temperature dependence of the two functions G(·) and R(·) is justified by experimental evidences demonstrating that m , KI and r values change with temperature (see [17]). Function m (Tp ) could be experimentally described by the following equation [18]: m (Tp ) = m,max T (Tp ),

Fig. 1. Schematic representation of the models implemented in this work. The arrows show the interconnection between the various models.

(5)

where m,max is the maximum value of m (Tp ) (s−1 ) and T (Tp ) is the temperature-dependent function reported in the following equation [18]:

R. De-Luca et al. / Journal of Process Control 55 (2017) 55–65 Table 1 Estimated values of the parameters used in the T function. Parameter

Description

Value

Tmin Tmax Topt r,max m,max KI,max

Minimum growth temperature Maximum growth temperature Optimum growth temperature Max. respiration coefficient Max. specific growth rate Max. half-saturation constant

−10.0 (◦ C) 42.1 (◦ C) 35.8 (◦ C) 2.01×10−6 (s−1 ) 6.48×10−5 (s−1 ) 7192.92 (W kg−1 )

T (Tp ) =

⎧ 0 ⎪ ⎪ ⎪ ⎨

57

as shown by [19]. A linear interpolation was carried out from these 6-h sampled data in order to get the meteorology at any time. The theoretical total amount of solar radiation reaching the external surface of the atmosphere at the location considered (H0 (t)) was derived from the computation of the solar angles by [20]: H0 (t) = Isc



1 + 0.033 cos

360N(t) 365



cos z (t)

if Tp ≤ Tmin (T − T

)(T − T

)

2

p max p min if Tmin < Tp ≤ Tmax ⎪ (Topt − Tmin )[(Topt − Tmin )(Tp − Tmin ) − (Topt − Tmax )(Topt + Tmin − 2Tp )] ⎪ ⎪ ⎩

0

= KI,max T (Tp ),

r (Tp )

= r,max T (Tp ).

(7)

Fitting these parameters was performed by using the Maximum Likelihood method included in the entity Parameter estimation of gPROMSTM software (4.1 version). The function T (Tp ) successfully fitted the experimental data for both m (Tp ), KI (Tp ) and r (Tp ). The parameter values are reported in Table 1. 2.2. Physical model of thermal evolution The temperature dynamics was obtained from a heat balance on the raceway pond (see [19]): dTp (t) = Qra,p (t) + Qra,s (t) + Qev (t) + Qconv (t) w V (t) cpw dt + Qcond (t) + Qi (t) + Qr (t) ,

(6)

if Tp > Tmax

The function T (Tp ) includes three parameters: Tmin , Tmax and Topt (◦ C). Tmin is the temperature below which the growth is assumed to be zero, Tmax is the temperature above which there is no growth nor respiration, Topt is the temperature at which m (Tp ) = m,max . Experimental values of m (Tp ), KI (Tp ) and r (Tp ) were extracted from the study of [17] conducted on Chlorella vulgaris as this model was shown to yield accurate predictions in outdoor photobioreactors. As r (Tp ) and KI (Tp ) exhibited similar evolution with temperature, the same function T (Tp ) was used for fitting the behavior of these two parameters at different temperatures: KI (Tp )

(9)

(8)

where cpw is the specific heat capacity of water (J kg−1 K−1 ), Qra,p (t) is the radiation flow from the raceway pond surface (W), Qra,s (t) is the total (direct+diffuse) solar irradiance (W), Qra,a (t) is the radiation flow from the air to the raceway pond system (W), Qev (t) is the evaporation heat flow (W), Qconv (t) is the convective heat flow at the raceway pond surface (W), Qcond (t) is the conductive heat flow with the ground at the raceway pond bottom (W), Qi (t) is the heat flow due to the water inflow (W), and Qr (t) is the heat flow associated with precipitation (W). The formal expression used to describe each heat flow in equation (8) was extracted from [19] and is recalled in Appendix 1. 2.3. Meteorological data Continuous weather data was extracted from the European Centre for Medium-Range Weather Forecast (ECMWF) website. This weather data was used to determine the dynamics of the air temperature Ta (t), the sky cloudiness CC(t), the relative humidity RH(t), the wind velocity vw (t) and the rain volumetric flux vr (t) as all these variables have a significant impact on raceway temperature

where Isc is the solar constant (1367 W m−2 ), N(t) is the integer variable representing the Julian day of the year and  z (t) is the zenith angle (the angle between a vertical axis at the location considered and the sun direction).  z (t) is given by: cos z (t) = sin  sin ı(t) + cos  cos ı(t) cos ω(t)

(10)

where ω(t) is the hour angle which varies linearly from - to over 24 h,  is the latitude and ı(t) is the solar declination, calculated as follows [20]; ı(t) = 23.45



2

284 + N sin 2

360 365



(11)

The clear-sky radiation Hc was given by [20]: Hc (t) = Hd,c (t) + HD,c (t) = (d,c (t) + D,c (t))H0 (t)

(12)

where  d,c (t) = Hd,c (t)/H0 (t) and  D,c (t)= HD,c (t)/H0 (t); Hd,c (t) and HD,c (t) are, respectively, the diffuse and the direct components of the clear-sky total irradiance (W m−2 ). The variable  D,c (t) and  d,c (t) were computed, respectively, through the correlations proposed in [21] and the Erbs correlation discussed in [20]. Then the solar irradiance at ground level Hs (t) was computed at any time taking into account the cloudiness CC(t), by using the Kasten and Czeplak correlation [21]:

0

Hs (t) =

if ω(t) < −ωs (t) or ω(t) > ωs (t) 3.4

Hc (t)(4 − 3(CC(t)/8) 4

)

(13) if − ωs (t) ≤ ω(t) ≤ ωs (t),

where CC(t) is the cloudiness value expressed in OKTAS (range 0−8), Hc (t) is the clear-sky total irradiance (W m−2 ). −ωs (t) and ωs (t) are, respectively, the hour angle values at sunrise and sunset calculated from the expression proposed by [20]. 2.4. System description The optimization strategy was applied on the weather data collected for the first 7.5 days of July 2012 in Nice, France (43◦ 42 11 N, 7◦ 15 57 E). The surface S of the open raceway was set at 100 m2 . The initial value of the raceway temperature Tp (t = 0) was set at the average of the air temperature Ta,avg over the considered week. Furthermore, the initial biomass concentration x was set at 0.4 kg/m3 and the initial raceway depth lp (t = 0) was set at 0.3 m. Finally, the inflow temperature Tin was set equal to the average air temperature Ta,avg . 3. MPC strategy 3.1. Objective criterion The control target was the maximization of microalgal productivity on a moving time window from the current time  n to  n + Th .

58

R. De-Luca et al. / Journal of Process Control 55 (2017) 55–65

Table 2 Computational cost evaluation for different discretisation times and horizons. The reference time Tr for a Samsung Electronics laptop, Intel(R) Core(TM) i7-3635QM CPU @ 2.40 GHz, 8.00 GB RAM, was Tr = 13.3 min. Control horizon

Control variables switching frequency 2h 1h 1/2 h

1 day

2 days

4 days

Tr 1.32 Tr 8.17 Tr

1.98 Tr 4.29 Tr 33.64 Tr

4.04 Tr 11.47 Tr 209.19 Tr

It consists in finding the optimal fresh medium volumetric flow rate (qin ) and the culture volumetric extraction rate (qout ). The control T

vector q(t) = (qin , qout ) was therefore the solution of the optimal control problem over the time horizon ( n ,  n + Th ):

 maxq(t) Pnet (t)

= maxq(t) s.t.

n +Th

(G( · ) − R( · ))V (t)dt

n

0 ≤ qi (t) ≤ qmax

,

(14)

lp min ≤ lp (t) ≤ lp max ˙ = g( , q(t), t), (0) = 0 where = (x, Tp )T is the state vector and g( , q(t), t) is the dynamical system coupling the biological and thermal dynamics (see Eqs. (1 and 8)). The upper bound of the two flow rates is qmax . The objective function Pnet (t) is the integral with respect to time of the difference between growth and respiration during the considered time window, therefore representing the net algal productivity. The computed optimal control is then applied for the period ( n ,  n+1 ). Since new weather forecasts are available at time  n+1 , a new control is then computed for the new period ( n+1 ,  n+1 + Th ). The procedure was iterated until the end of the cultivation period (at day 7.5). We assume that temperature and biomass density are perfectly on-line measured and known at the beginning of each new optimization period (biomass can for example be assessed by simple measurements of the optical density [22]). The study in [13] addressed the same optimization problem but imposed the pond depth to be constant, which significantly constrained the evolution of qin (t) and qout (t). In the present study, the pond depth was constrained within a minimal (0.05 m) and a maximal (0.5 m) value, which turns out to have a significant impact on system temperature as discussed in section 4. The optimization task was implemented through gPROMSTM (4.1 version) default solver NLPSQP, which uses a sequential quadratic programming (SQP) method to solve nonlinear programming (NLP) problems (optimization tolerance is set to 0.001). The two control inputs qin (t) and qout (t) were numerically implemented as piecewise constants within the range [0–1] m 3 s−1 . Each new optimization problem is initialized with the flow vector resulting from the previous optimization which both reduced the computation time and decreases the risk of local minima. Table 2 shows the computational times required for different optimization tasks. The reported results were obtained through the variation of the switching frequency of each input variable (30 min, 1 h, 2 h) and the implementation of different control horizons (1 day, 2 days and 4 days). We selected a reasonable trade-off between accuracy and computation time with a 1 h input discretization time for a 4 day horizon. This required 2.5 h on our computer, which is realistic for an on-line implementation. Each case study was initialized by setting the two control variables (qin (t) and qout (t)) equal to 10−4 m3 s−1 .

Fig. 2. Global and local optimal profiles of: (a) qin (t) and (b) qout (t). Global optimization (red line) represents the profiles obtained by taking into account the time horizon ( 1 ,  f ) whereas the local optimal profiles (black lines) represent the results obtained for all the sub-optimization tasks conducted in the intervals ( n ,  n + Th ), with n=0:6 and Th =4. (c) shows the profiles obtained with the global optimization (solid line) and the standard control of open ponds [23,24] (dashed line). The background is colored in white at daytime and in grey at nighttime. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

R. De-Luca et al. / Journal of Process Control 55 (2017) 55–65

59

Fig. 3. Three-days zoom of (a) qin and (b) qout optimized profiles. (The background is colored in white at daytime and in grey at nighttime.)

4. Optimization with perfect weather forecasts

4.2. Interpretation of the optimal control strategy

4.1. Comparison with a reference management strategy

Fig. 2a and b present qin (t) and qout (t) profiles for the entire duration of the cultivation. The two graphics show a qualitatively daily periodic pattern for both input variables, despite the weather variability between different days of cultivation. The only exception was on the first day of cultivation at nighttime, when a peak of qout (t) was triggered in order to rapidly minimize the initial high biomass content in the system. As a consequence, and for the sake of clarity, the majority of the results will be reported only for three days of cultivation (between t = 2 days and t = 5 days). Fig. 3a and b detail the optimal inflow rate and the outflow rate profiles for three days of cultivation. The following discussion focuses on the analysis of the optimal control strategy and investigates how it significantly contributes to increase net productivity. Fig. 5a–c show, respectively, the pond depth lp profile, the pond temperature Tp dynamics and the biomass concentration x. The analysis of the results is split in four phases, from morning to night:

This first study was based on a perfect weather forecast. In Fig. 2a and b the optimal variable profiles obtained by considering both the global optimization ( 1 ,  f ) and the recurrent sequence of suboptimization tasks defined in section 3 (with Th set equal to 4 days) are reported. No significant variations between the two approaches were registered, except for the final peaks of qin (t) and qout (t) at the end of the final day of each interval ( n ,  n + Th ). Such behavior is related to the fact that the 4th day of each time window is considered by the optimizer as the last day of cultivation. A complete emptying of the open raceway pond was therefore required to collect all the residual biomass in the system. However, since the optimization update was carried out every day, no practical effect of such peaks on the raceway operation was registered and almost a perfect superimposition of the global and local optimization curves was obtained. The net productivity Pnet (t) obtained under optimal operation is reported in Fig. 2c. It shows that the cumulative productivity for 7.5 days in Nice was 28.2 kg, which, when extrapolated, corresponds to a productivity of 137.2 t ha−1 year−1 . This high value was mostly due to the fact that Nice benefits from a sunny climate and temperature values close to Topt at daytime in summer (see Fig. 5b). The productivity at this period of the year is therefore close to the maximum achievable productivity in Nice and would be much lower during winter. The assessment of the gain of productivity was carried out through a reference simulation under ‘standard conditions’. These standard conditions correspond to a typical management of microalgae outdoor ponds: constant depth of 0.25 m and dilution rate of 0.1 day−1 ([23], [24]) for the same climatic conditions. In these standard conditions, the productivity during the period was 13.2 kg. This highlights a dramatic drop in productivity when no control is carried out during hot days, leading to culture overwarming and high algal concentrations. This factor 2.13 in the final productivity compared to standard operation illustrates the striking benefit of accounting for weather forecast to determine the raceway pond inflow and outflow. Ideally, the optimal control problem should be solved on the interval ( n , Tf ). However, the dimension of this problem when Tf −  n is larger than 5 days induces a high computational cost, and the risk of non convergence due to non convexity. Of course, the larger Th the better, but the choice of Th is practically determined by the computation time to solve the problem on a larger interval. It turns out that (see Fig. 2a and b) Th = 4 days leads to an efficient trade-off (see Table 2), while considering a larger time window did hardly influence the optimal solution on the interval ( n ,  n+1 ).

• Morning: Fig. 5(a) shows that no water is injected into or extracted from the open pond in the morning. The pond depth was therefore kept constant and at a relatively low value (see Fig. 5(a)). Reducing the raceway pond depth in the morning decreases the thermal inertia, and thus accelerates the temperature increase in order to rapidly reach the optimal temperature between 9 and 12 am (see Fig. 5(b)). Interestingly, the pond depth was slightly higher than the minimal allowed pond depth (0.05 m) simply because removing more culture from the open pond would lead to further decrease the biomass content in the pond, hence the productivity. • Afternoon: In this phase qin (t) followed a ‘bell curve’ profile from about noon until late afternoon, whereas qout (t) followed the same dynamics with a slight delay (see Fig. 3(a) and (b)). Fig. 5(b) shows that the pond temperature Tp (t) during the afternoon was kept near its optimal value Topt (35.8 ◦ C) in spite of high solar flux and high air temperatures in the middle of the afternoon. The good level of temperature control can therefore be explained by the additional action of replacement of the raceway pond water by colder fresh medium. Furthermore, Fig. 5(a) shows that the raceway depth lp (t) increased until mid-afternoon and then decreased. This dynamics was consistent with the delay of the bell profile of qout (t) with respect to qin (t) profile (Fig. 3(a) and (b)). The consequent increase of depth enhanced the thermal inertia of the system, therefore limiting temperature deviation from Topt . The depth increase during the afternoon was an additional strategy used by the optimizer to maintain the raceway temperature at its optimal value when the ‘flushing’ strategy discussed above was not satisfactory in terms of temperature control.

60

R. De-Luca et al. / Journal of Process Control 55 (2017) 55–65

For example, if the function fcomp (t) is higher than 1, it means that net productivity could be improved by increasing the amount of biomass in the system. Conversely, values lower than 1 indicate that the respiration rate is higher than the growth rate. Consequently, diluting the system would increase productivity. The optimal biomass concentration xopt (t) is therefore reached by keeping fcomp (t) equal to 1. Fig. 4 illustrates the fcomp (t) profile for the same three days of cultivation discussed before. • Sunset: Fig. 3(a) and (b) show that a fraction of the culture was extracted at sunset without any replacement by fresh medium. The optimal strategy in this case consisted of decreasing the raceway pond depth (Fig. 5(a)). In addition to biomass removal, decreasing the depth indeed allowed a fast decrease of temperature, hence lowering the respiration rate. • Night: At nighttime the depth was maintained at its sunset value all night long (see Fig. 5(a)). 5. Accounting for uncertain weather forecasts Fig. 4. Three-days zoom of fcomp (t) optimized profile. (The background is colored in white at daytime and in grey at nighttime.)

On top of this temperature control strategy, the optimizer actions led to maintain the algal concentration at its optimal value during the afternoon. Indeed, it has been shown that, for constant [25] or varying [13] light, the ratio between the growth and the respiration rates at the raceway bottom should be close to one to ensure optimal productivity. This ratio can be computed using the so-called compensation function fcomp (t), expressed through the following equation: fcomp (t) =

m,max (H Hs (t)e−x(t)lp (t) /KI,max T (t) + H Hs (t)e−x(t)lp (t) ) . (15) r,max

5.1. Numerical implementation of weather forecast uncertainty Despite the increase in meteorological accuracy and computing power, the weather prediction capability is in the range of a few days. To simulate weather forecast uncertainty, we assumed a perfect forecast during the first day and we simply shifted the meteorological data by one day for the rest of the period. This is probably a strong perturbation, in comparison with the accuracy of the existing meteorological models (see [26,27] for a detailed study) but it generates a credible weather forecast. This is illustrated on Fig. 6 for cloudiness data, but the logic was used for all the five weather variables extracted by ECMWF website.

Fig. 5. Three-days zoom of (a) lp (t), (b) Tp (t) and (c) x(t) for the trajectory resulting from the optimal control. (The background is colored in white at daytime and in grey at nighttime.)

R. De-Luca et al. / Journal of Process Control 55 (2017) 55–65

Fig. 6. The real and predicted profiles for cloudiness are reported. The solid line represents the real weather dynamics whereas the dashed line represents the weather forecast used for the optimization. (The vertical dashed-dotted line represents the artificial division between unaltered and modified weather data for a generic ( n ,  n + Th ) time window.)

More accurately, for each preset optimization time window ( n ,  n + Th ) (where  n is the current time and Th is the length of the considered time horizon) the weather forecast data was corrupted in the following way: • No discrepancies between real and predicted data were introduced in the time window ( n ,  n+1 );

61

Fig. 8. Five-days zoom of the CC profiles. The solid black line represents the real weather data, whereas the black and red dashed lines represent, respectively, the weather forecast for the first and the second optimization windows. The horizontal arrows and the dashed-dotted vertical lines are used to point out each optimization window. The first window ( 1 ,  1 + Th ) is in black and the second ( 2 ,  2 + Th ) is in red (the background is colored in white at daytime and in grey at nighttime.) (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

• For each jth day in the interval ( n+1 ,  n + Th ), the data used for the optimization were the same collected for the (j + 1)th day; as an example, if we consider the first 4 days of cultivation, the data used for day 2 were the same as the real ones collected for day 3, and so on till day 4.

Fig. 7. (a) and (b) show the optimal control values of qin (t) and qout (t) by taking into account the uncertain weather forecast data reported in Fig. 6, whereas in (c) and (d) the effect of forecast uncertainty on Tp and net productivity is represented. (The background is colored in white at daytime and in grey at nighttime. The arrows highlight the negative effects of wrong forecasts on net productivity as a consequence of qin (t) and qout (t) optimal profiles.)

62

R. De-Luca et al. / Journal of Process Control 55 (2017) 55–65

5.2. Impact of inaccurate weather on the optimal solution over the ( n ,  n + Th ) time window In this paragraph we assess the impact of inaccurate weather forecasts. Optimizing the control strategy by setting large time windows and, consequently, using highly uncertain weather forecasts, has a great impact on the net productivity of the system. Fig. 7(d) shows that at day 3 there is a loss of productivity with respect to the predicted profile. This critical point affected the productivity of the 4-day time window by lowering its value of 3.2 kg with respect to the ideal case previously described. This result can be explained by considering the qin (t) and qout (t) optimal profiles. Focusing on day 3, Fig. 7(a) and (b) shows low control peaks due to the fact that a cloudy day was wrongly predicted. Since day 3 was a sunny day the optimal control strategy was not sufficient to cope with high temperatures and high solar irradiances. In fact, Fig. 7(c) shows that at day 3 the raceway pond temperature overcame the Tmax value hindering the growth rate of the culture. In the next paragraph, we examine how a daily update of the weather forecasts together with a 4 day time window can compensate for the possible impact of a wrong meteorology. 5.3. MPC accounting for inaccurate meteorological forecasts In this section, the predictive controller is implemented with Th = 4 days by using the MPC approach described in Section 3.

Fig. 10. (a) lp (t), (b) Tp (t) and (c) x(t) optimal control values. The horizontal arrows and the dashed-dotted vertical lines are used to point out each optimization window. For the sake of simplicity, only the first two optimization windows ( 1 ,  1 + Th ) (black) and ( 2 ,  2 + Th ) (red) are represented (The background is colored in white at daytime and in grey at nighttime.) (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

Fig. 9. (a) qin (t) and (b) qout (t) optimal control values. The horizontal arrows and the dashed-dotted vertical lines are used to point out each optimization window. For the sake of simplicity, only the first two optimization windows ( 1 ,  1 + Th ) (in black) and ( 2 ,  2 + Th ) (in red) are represented (The background is colored in white at daytime and in grey at nighttime.) (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

For the sake of simplicity the results are reported only for the first two optimization iterations. The inflow rate qin (t) and the outflow rate qout (t) profiles are reported in Fig. 9(a) and (b) whereas the raceway pond depth lp (t), the raceway pond temperature Tp (t) and the biomass concentration x(t) are reported in Fig. 10. The MPC approach with daily re-optimization allows to maintain productivity values close to the ideal case (28 kg). These results

R. De-Luca et al. / Journal of Process Control 55 (2017) 55–65

63

demonstrate that the optimal control strategy is mainly linked to weather conditions in the first days of optimization and that updating the weather forecast only once per day is sufficient to prevent the inaccuracies caused by weather forecast uncertainty.

where Lw is the water latent heat (J kg−1 ) and me (t) is the mass flux caused by evaporation (kg m−2 s−1 ), given by the equation:

6. Conclusions

where K is the mass transfer coefficient (m s−1 ), Pw (t) and Pa (t) are, respectively, the saturated vapor pressure (in Pa) at Tp (t) and Ta (t), calculated through the empirical correlation:

The main result of our study is the key role played by temperature control for optimizing microalgal productivity. Several studies have focused on control strategies deriving an optimal microalgal density, balancing growth and respiration at the bottom of the reactor (compensation condition). The MPC strategy however gives to this target only the second rank. It first adapts the raceway pond volume and flow rates in order to get as close as possible to the optimal temperature for growth, and always prevents the temperature from exceeding this optimal value. The analysis of the rational behind the optimal control revealed that this productivity gain was achieved via two main mechanisms: ‘flushing’ the culture and controlling the raceway pond depth. The rather periodical control actions give the pace for inflow and outflow. The compensation condition is then tracked once temperature is close to the optimal value. This control scheme turns out to significantly increase the productivity of algal open raceway ponds compared to standard operation with constant depth and dilution rates. The proposed control strategy is highly dependent on the weather forecasts, and an accurate prediction must be available at least for one day to guarantee an efficient process optimization. Inaccurate predictions can be compensated by a more frequent optimization update. It is shown in [28] that the complex optimal strategy can be reduced to six simple rules whose implementations are simpler. Using these simple rules to optimize raceway operation has therefore the potential to i. minimize computational costs, ii. supply guidelines for practical operation of open raceway ponds and iii. to develop simple closed loop control loops in order to reduce the sensitivity of the control strategy to inaccurate weather forecasts. Acknowledgements

Appendix A. Details of the thermal model Assuming the water surface to be gray-diffuse, the radiation Qra,p (t) from the raceway pond surface to the atmosphere was given by Stefan–Boltzmann’s law: (A.1)

where εw is the emissivity of water and  SB is the Stefan–Boltzmann constant (W m−2 K−4 ). The heat flow associated with solar irradiance Qra,s (t) was described as: Qra,s (t) = (1 − fa )Hs (t)S,

(A.2)

where fa is the fraction of the photosynthetically active radiation (PAR) that is effectively used during the photosynthetic process. The heat flow generated by air radiation was expressed through the equation: Qra,a (t) = εa εw SB Ta4 (t) S,

(A.3)

where εa is the emissivity of the air and Ta (t) is the air temperature (K). The evaporation heat flow Qev (t) was given by the equation: Qev (t) = −me (t)Lw S,

me (t) = K

Pw (t) RH · Pa (t) − Tp (t) Ta (t)



Pi (t) = 3385.5e



Mw , R

(A.5)

−8.0929+0.97608(Ti (t)+42.607−273.15)0.5

,

(A.6)

where the i index represents air or water, RH is the relative air humidity above the raceway pond surface, Mw is water molecular weight (kg mol−1 ) and R is the universal ideal gas constant (Pa m3 mol−1 K−1 ). The mass transfer coefficient K in (A.5) was evaluated through the following correlation:



Sh =

0.035Re0.8 Sc 1/3 0.628Re

0.5

Sc

1/3

for turbulent flows (A.7) for laminar flows

where Sh = KL/Dw,a , Re = Lvw (t)/a and Sc = a /Dw,a . Sh, Re and Sc are, respectively, the Sherwood, Reynolds and Schmidt dimensionless numbers, L is the characteristic raceway length (m), Dw,a is the mass diffusion coefficient of water vapor in air (m2 s−1 ), vw (t) is the wind velocity (m s−1 ) and a is the kinematic viscosity of the air (m2 s−1 ). The convective flow Qconv (t), defined by the equation: Qconv (t) = hconv (Ta (t) − Tp (t))S,

(A.8)

was calculated by evaluating the heat transfer coefficient hconv (W m−2 K−1 ) through the following correlation:



Nu =

0.035Re0.8 Pr 1/3

for turbulent flows,

0.628Re0.5 Pr 1/3

for laminar flows,

(A.9)

where Nu =hconv L/a and Pr =a /˛a . Nu and Pr are, respectively, the Nusselt and the Prandtl dimensionless numbers, a is the air thermal conductivity (W m−1 K−1 ) and ˛a the air thermal diffusivity Table A1.1 Parameter values of the thermal model

The authors acknowledge the support of the ANR projects Purple Sun (ANR-13-BIME-004) and Phycover (ANR-14-CE04-0011). RDL and FB acknowledge the support of project PRIN 2012 (2012XSAWYM).

Qra,p (t) = −εw SB Tp4 (t) S,



(A.4)

Water parameters water density w water heat capacity c pw water latent heat Lw εw water emissivity water molecular weight Mw

998 (kg m−3 ) 4.18×103 (Jkg−1 K−1 ) 2.45106 (J kg−1 ) 0.97 (–) 0.018 (kg mol−1 )

Soil parameters soil thermal conductivity ks soil heat capacity c ps soil density s soil temperature at lsref = 4.5m Tsref

1.7 (W m−1 K−1 ) 1.25×103 (J kg−1 K−1 ) 1.9×103 (kg m−3 ) 286.75 (K)

Air parameters air emissivity εa air kinematics viscosity a air thermal conductivity a air thermal diffusivity ˛a mass diffusion coefficient of Dw,a water vapor in air

0.8 (–) 1.5×10−5 (m s−1 ) 2.6×10−2 (W m−1 K−1 ) 2.2×10−5 (m s−1 ) 2.4×10−5 (m s−2 )

Raceway pond parameters initial raceway pond volume V0 raceway surface S raceway characteristic length L power law exponent ˛ wind velocity height z wind sensor height z0 algal absorption fraction fa

30 (m3 ) 100 (m2 ) 10 (m) 0.29 (–) 0.5 (m) 10 (m) 2.5 (%)

Universal constants Stephan–Boltzmann constant  SB ideal gas constant R

5.67×10−8 (W m−2 K−4 ) 8.314 (Pa m3 mol−1 K−1 )

64

R. De-Luca et al. / Journal of Process Control 55 (2017) 55–65

Fig. A2.1. Meteorological data for July 2012 in Nice. The background is colored in grey at daytime and in white at nighttime).

R. De-Luca et al. / Journal of Process Control 55 (2017) 55–65

(m2 s−1 ). Since meteorological stations measure the wind velocity at a height z0 (m) which usually differs from the height at which vw (t) is needed in the previous correlations (0.5 m), the conversion of the wind velocity from the height z0 to 0.5 m was performed by using the following expression:

 z ˛

vw (t) = v0 (t)

[6]

(A.10)

[7]

where v0 (t) is the wind velocity (m s−1 ) measured at height z0 and ˛ is a power law exponent. The conductive heat flow between the raceway and the soil was described through the Fourier’s law-based equation:

[8]

z0

Qcond (t) = ks S

,

[5]

dTs (t, z) (z = 0) , dz

(A.11)

where ks is the soil conductivity (W m−1 K−1 ) and Ts (t) is the soil temperature (K). The soil temperature Ts (t, z) was obtained from the following equation and boundary/initial conditions: cps s

dTs (t, z) dt

[9] [10]

[11]

[12]

2

=

ks

d Ts (t, z) dz 2

[13]

s.t. Ts (t, z = 0)



Ts t, z = lsref d2 Ts dz

2

= Tp (t)

(t = 0, z)

=

[14]

,

(A.12)

Tsref

[15]

[16]

= 0

[17]

where cps is specific heat capacity of the soil (J kg−1 K−1 ), s is the soil density (kg m−3 ), and Tsref is the soil temperature (K) at the reference depth lsref (m) reported in Table A1.1. The heat flow associated with fresh medium inflow Qin (t) was modeled using the equation: Qi (t) = w cpw qin (t)(T in − Tp (t)),

[19]

(A.13)

where Tin is the water inflow temperature (K). Finally, the rain heat flow Qr (t) was expressed as: Qr (t) = w cpw vr (t)(Ta (t) − Tp (t))S.

[18]

[20] [21]

(A.14)

All the parameters values used in the energy balance (extracted from [19]) are reported in Table A1.1.

[22]

[23]

Appendix B. Meteorological data The real weather data extracted from the European Centre for Medium-Range Weather Forecast (ECMWF) website were: the air temperature Ta (t), the sky cloudiness CC(t), the relative humidity RH(t), the wind velocityvw (t) and the rain volumetric flux vr (t) during the simulated cultivation period. Fig. A2.1 represents the profiles of the implemented meteorological data. References [1] T.M. Mata, A.A. Martins, N.S. Caetano, Microalgae for biodiesel production and other applications: a review, Renew. Sustain. Energy Rev. 14 (1) (2010) 217–232, http://dx.doi.org/10.1016/j.rser.2009.07.020. [2] K. Skjånes, C. Rebours, P. Lindblad, Potential for green microalgae to produce hydrogen, pharmaceuticals and other high value products in a combined process, Crit. Rev. Biotechnol. 33 (2) (2013) 172–215. [3] P.M. Foley, E.S. Beach, J.B. Zimmerman, Algae as a source of renewable chemicals: opportunities and challenges, Green Chem. 13 (6) (2011) 1399–1405. [4] F. Grognard, A.K. Akhmetzhanov, O. Bernard, Optimal strategies for biomass productivity maximization in a photobioreactor using natural light,

[24]

[25]

[26]

[27]

[28]

65

Automatica 50 (2014) 359–368, http://dx.doi.org/10.1016/j.automatica.2013. 11.014. E.F. Camacho, C.B. Alba, Model Predictive Control, Springer Science & Business Media, 2013. F. Oldewurtel, A. Parisio, C.N. Jones, D. Gyalistras, M. Gwerder, V. Stauch, B. Lehmann, M. Morari, Use of model predictive control and weather forecasts for energy efficient building climate control, Energy Build. 45 (2012) 15–27. Y. Zong, D. Kullmann, A. Thavlov, O. Gehrke, H.W. Bindner, Application of model predictive control for active load management in a distributed power system with high wind penetration, IEEE Trans. Smart Grid 3 (2) (2012) 1055–1062. M. Berenguel, F. Rodrıguez, F. Acién, J. Garcıa, Model predictive control of ph in tubular photobioreactors, J. Process Control 14 (4) (2004) 377–387. S. Tebbani, F. Lopes, R. Filali, D. Dumur, D. Pareau, CO2 Biofixation Bymicroalgae: Modeling, Estimation and Control, Wiley ISTE, 2014. S. Tebbani, F. Lopes, R. Filali, D. Dumur, D. Pareau, Nonlinear predictive control for maximization of CO2 bio-fixation by microalgae in a photobioreactor, Bioprocess Biosyst. Eng. 37 (1) (2014) 83–97. O. Bernard, Hurdles and challenges for modelling and control of microalgae for CO2 mitigation and biofuel production, J. Process Control 21 (2011) 1378–1389, http://dx.doi.org/10.1016/j.jprocont.2011.07.012. Q. Béchet, A. Shilton, B. Guieysse, Modeling the effects of light and temperature on algae growth: state of the art and critical assessment for productivity prediction during outdoor cultivation, Biotechnol. Adv. 31 (8) (2013) 1648–1663, http://dx.doi.org/10.1016/j.biotechadv.2013.08.014. ˜ R. Munoz-Tamayo, F. Mairet, O. Bernard, Optimizing microalgal production in raceways systems, Biotechnol. Prog. 29 (2) (2013) 543–552, http://dx.doi.org/ 10.1002/btpr.1699. O. Bernard, A.C. Boulanger, M.O. Bristeau, J. Saint-Marie, A 2D model for hydrodynamics and biology coupling applied to algae growth simulations, ESAIM: Math. Numer. Anal. 47 (5) (2013) 1387–1412. J. Mendoza, M. Granados, I.D. Godos, F. Acién, E. Molina, C. Banks, S. Heaven, Fluid-dynamic characterization of real-scale raceway reactors for microalgae production, Biomass Bioenergy 54 (2013) 267–275. F.J. Weissing, J. Huisman, Growth and competition in a light gradient, J. Theor. Biol. 168 (3) (1994) 323–336. Q. Béchet, P. Chambonnière, A. Shilton, G. Guizard, Algal productivity modeling: a step toward accurate assessments of full-scale algal cultivation, Biotechnol. Bioeng. 112 (5) (2015) 987–996, http://dx.doi.org/10.1002/bit. 25517. O. Bernard, B. Remond, Validation of a simple model accounting for light and temperature effect on microalgae growth, Bioresour. Technol. 123 (2012) 520–527, http://dx.doi.org/10.1016/j.biortech.2012.07.022. Q. Béchet, A. Shilton, J.B.K. Park, R.J. Craggs, B. Guieysse, Universal temperature model for shallow algal ponds provides improved accuracy, Environ. Sci. Technol. 45 (2011) 3702–3709, http://dx.doi.org/10.1021/es1040706. J.A. Duffie, W.A. Beckman, Solar Engineering of Thermal Processes, Little, Brown & Co., Boston, 1958. T.R. Marthews, Y. Malhi, H. Iwata, Calculating downward longwave radiation under clear and cloudy conditions over a tropical lowland forest site: an evaluation of model schemes for hourly data, Theor. Appl. Climatol. 107 (3-4) (2012) 461–477, http://dx.doi.org/10.1007/s00704-011-0486-9. I. Havlik, P. Lindler, T. Sheper, K.F. Reardon, Online monitoring of large cultivations of microalgae and cyanobacteria, Trends Biotechnol. 31 (7) (2013) 406–414, http://dx.doi.org/10.1016/j.tibtech.2013.04.0. O. Jorquera, A. Kiperstok, E.A. Sales, M. Embiruc¸u, M.L. Ghirardi, Comparative energy life-cycle analyses of microalgal biomass production in open ponds and photobioreactor, Bioresour. Technol. 101 (2010) 1406–1413, http://dx. doi.org/10.1016/j.biortech.2009.09.038. J.N. Rogers, J.N. Rosenberg, B.J. Guzman, V.H. Oh, L.E. Mimbela, A. Ghassemi, M.J. Betenbaugh, G.A. Oyler, M.D. Donohue, A critical analysis of paddlewheel-driven raceway ponds for algal biofuel production at commercial scales, Algal Res. 4 (2014) 76–88, http://dx.doi.org/10.1016/j. algal.2013.11.007. H. Takache, G. Christophe, J.F. Cornet, J. Pruvost, Experimental and theoretical assessment of maximum productivities for the microalgae Chlamydomonas reinhardtii in two different geometries of photobioreactors, Biotechnol. Prog. 26 (2010) 431–440, http://dx.doi.org/10.1002/btpr.356. E. Lorenz, J. Hurka, D. Heinemann, H.G. Beyer, Irradiance forecasting for the power prediction of grid-connected photovoltaic systems, IEEE J. Sel. Top. Appl. Earth Observ. Remote Sens. 2 (1) (2009) 1939-1404, http://dx.doi.org/ 10.1109/JSTARS.2009.2020300. R. Perez, E. Lorenz, S. Pelland, M. Beauharnois, G. Van-Knowe, D. Heinemann, J. Remund, S.C. Müller, W. Traunmüller, G. Steinmauer, D. Pozo, J.A. Ruiz-Arias, V. Lara-Fenego, L. Ramirez-Santigosa, M. Gaston-Romero, L.M. Pomares, Comparison of numerical weather prediction solar irradiance forecasts in the US, Canada and Europe, Sol. Energy 94 (2013) 305–326, http://dx.doi.org/10.1016/j.solener.2013.05.005. R. De-Luca, Q. Béchet, F. Bezzo, O. Bernard, Meteorological data-based optimal control strategy for microalgae cultivation in open pond systems, AIChE Journal. (submitted for publication).