Weather forecast error modelling and performance analysis of automatic greenhouse climate control

.


Introduction
The Netherlands is one of the biggest exporters of vegetables in the world, but this comes at a price.In 2018, the Dutch horticultural industry consumed 100:5 PJ, of which only 7:4 PJ were produced in a sustainable manner.This use of energy resulted in a CO 2 emission of 5:7 Mt (Velden & Smit, 2019).The Dutch horticultural industry signed an agreement with the Dutch government to decrease the CO 2 emission and its environmental footprint.Among other innovations, automatic control of the greenhouse is likely to contribute to achieving the goals set in the Dutch agreement and may lead to a more sustainable cultivation worldwide (van Straten & van Henten, 2010) heat stored in heat buffer, J:m À2 x t the state of the system at the present time Automatic control can be used to compute suitable trajectories for the control inputs to the greenhouse system (e.g. level of operation of the combined heat and power unit (CHP), ventilation rate, CO 2 injection) in order to improve, for example, yield, CO 2 emission or a combination of both.To realise optimal control, the literature often proposes optimal control (Kuijpers et al., 2021;Seginer et al., 2018;Tap, 2000;van Straten & van Henten, 2010) among various other control approaches.In optimal control, the control action is the result of an objective function optimisation over a future time interval.Throughout this time interval, the operation of the greenhouse system is affected by various weather variables, such as global radiation, temperature and wind speed.Optimisation over a future time interval thus requires a forecast of the relevant weather variables.Due to necessary approximations in numerical weather models used for weather forecasting and the chaotic nature of the weather system, weather forecasts may not match the actual realisation of the weather variable, introducing uncertainty in the weather forecast.The uncertainty of a weather forecast typically increases with forecast lead time.The forecast lead time is the period of time between the instance of time when the forecast is published and the time instance at which the weather variable is predicted.
In the literature on optimisation-based greenhouse climate control t, various approaches to modelling and incorporating weather forecasts are presented.When simulating past time instances, the actual realisation of the weather can be substituted for the weather forecasts (Achour et al., 2020;Ferreira & Ruano, 2008;Seginer et al., 2017b;van Beveren et al., 2015).Weather forecasts can also be synthesised by, e.g.assuming periodic weather (Seginer et al., 2017a) combined with a stochastic variable (Chen et al., 2018) or by assuming constant weather (Ioslovich & Seginer, 2002).Irrespective of the approach used, most published research assumes that the weather forecast matches the actual realisation of the weather variable, i.e. no weather forecast error.The lazy-man weather forecast presented in Tap (2000) and used in van Ooteghem (2007), however, does include a weather forecast error.Also, in documented experimental studies, such as Ghoumari et al. (2005), a weather forecast error was included through the use of weather forecasting services.In Su et al. (2021), the weather forecasts are based on data collected in previous years.The latter researches, alongside most other documented researches in which a weather forecast error is present, did not incorporate the uncertainty of the weather forecast error in the synthesis of the controller.Chen and You (2021) presents an approach in which the uncertainty in the weather forecasts is included, this approach is, however, limited to short prediction horizon.The latter makes it impossible to use an economic cost function as presented in the precursor of the work presented here (Kuijpers et al., 2021).
If the uncertainty in the weather forecast error is not included in the synthesis of the controller, the resulting control action may be optimal with respect to the employed objective function for some realisations of the weather but suboptimal for others, i.e. risk is not taken into account (Doeswijk, 2007).In turn, a system takes risk into account if it satisfies robust performance conditions, i.e. performance specifications are met for a predetermined set of possible realisations of the uncertainty.In Mayne (2014) and Saltık et al. (2018), various model predictive control (MPC) algorithms are discussed that are tailored to uncertain systems.If the weather forecast error is considered as a stochastic variable, various algorithms can be employed to obtain robust performance.The research presented in this paper builds upon the work presented in Oldewurtel et al. (2014) and Zhang et al. (2013), in which stochastic MPC and randomised MPC (RMPC) algorithms were considered for building climate control.The latter control algorithms may also be viable solutions to the challenges in the greenhouse control problem, primarily as we hypothesise that the uncertainty in weather forecasts considerably affects the performance of both applications.
To the best of our knowledge it is unclear to what extent the weather forecast error affects the performance of an optimally controlled greenhouse system.This research aims to quantify the loss in performance of the system due to weather forecast errors when the controller synthesis neglects the effect of the uncertainty.The performance of the greenhouse system is measured in terms of the operational return, which balances the cost of resources (resource use Â cost) with the income through yield (yield Â product price).This research employs weather forecasts that are synthesised using a stochastic weather forecast error model, inspired by Oldewurtel et al. (2014) and Zhang et al. (2013).The forecast error model is based on an autoregressive (AR) model with a stochastic element.Here, a model is used to predict the weather forecast error, in contrast to Huang and Chalabi (1995), who showed an approach in which the windspeed is predicted using an adaptive AR model.The parameters of the forecast error model are identified using historical observations and historical forecasts from the X set of admissible values for the states y c the interaction between the crop and the greenhouse y g the interaction between the greenhouse and the crop KNMI, more specifically, forecasts from the HARMONIE-AROME Cy40 model (Bengtsson et al., 2017).Similar to the KNMI, the synthesised weather forecasts are published every 6 h.A Kalman filter updates the previously published forecast until a new forecast is published, using the measured weather at the greenhouse site, an approach similar to Doeswijk (2007) and Oldewurtel et al. (2014).This research also aims to quantify the loss in performance of the system due to weather forecast errors when the controller synthesis takes into account the effect of the uncertainty.To this extent, the RMPC algorithm presented in Zhang et al. (2013) was also evaluated.The contribution of this research is three-fold:

Acronyms
The parameters in the stochastic error model proposed in Oldewurtel et al. (2014) and Zhang et al. (2013) are identified based on historical data from the HARMONIE-AROME Cy40 model for the Dutch climate.
The loss in performance of the optimally controlled greenhouse system due to weather forecasts is quantified for non-zero weather forecast error when the controller neglects the effect of the uncertainty.
The loss in performance of the optimally controlled greenhouse system due to weather forecasts is quantified when controlled with RMPC, an algorithm in which the effect of the uncertainty is explicitly included.
The paper is the first to investigate the impact of weather forecast errors on the optimal control of greenhouse which so far has been tackled without taking the uncertainty in these forecast errors into account.We conclude that forecast error have a significant impact on control inputs but not on return.These findings provide useful insights and play an important role in the area of optimal greenhouse control.
The remainder of this paper is organised as follows, Section 2 elaborates on the models employed for the greenhouse and the weather forecasts, as well as the various control approaches that are used.The simulation studies and corresponding results are presented in Section 3 and the results are discussed in Section 4. Directions for future work and the conclusion of this research are presented in Section 5.

Models & methods
The research presented here builds upon the greenhouse control problem as presented in Kuijpers et al. (2021)

Greenhouse control problem
The controlled greenhouse system is graphically represented by the block diagram in Fig. 1.The model of the greenhouse system is composed of the energy management system model S E , greenhouse climate system model and lighting system model S G and crop growth and transpiration modelS C .The greenhouse system is modelled with a state-space representation with states x2R nx , controllable inputs u2R nu and uncontrollable inputs d2R n d .The components of the state vector x and controllable inputs vector u are listed in Tables 1 and 2, respectively.The interaction between the greenhouse climate model and the crop model (temperature, CO 2 concentration, radiation and relative humidity) and vice versa (assimilation and transpiration) are denoted by y g and y c in Fig. 1, respectively.The dynamical model The continuous-time controllable inputs to the greenhouse system u can be updated every 15 min, t l ¼ 15,60 ¼ 900 s and are held constant in between samples, following Kuijpers et al. (2021).These control inputs are typically inputs to other control systems realising a desired setpoint, e.g.e.g. the heating system is controlled using u boi 2R, in practice, instead of a control, this will be a setpoint provided to a low-level control system.For the sake of computational efficiency, these inputs are discretised using a zero-order hold with sampling time t l .In the simulation studies presented here, the discrete-time control inputs to the greenhouse system u d 2R nu resulted from an optimisation problem solved by the receding horizon optimal controllerS M , the resulting control inputs are denoted by u * d 2R nu .The controller aimed to optimise the operational return J2R (V m À2 ),  2. These constraints represent either equipment limits or limits to the domains of the models being used.The authors avoid using soft constraints as these require the tuning of artificial parameters and will dilute the economic cost function.All inputs have a continuous domain.The n e inequality constraints are expressed in terms of functions h : R nx Â R nu /R ne , these are lower and upper bounded by q l 2R ne and q u 2R ne , respectively.The optimisation problem aims to find u * d by maximizing the operational return within the feasible region outlined by the constraints The state of the system at the present time is represented by x t 2R nx .Vector x d ðjjkÞ is the predicted1 state at future time Fig. 1 e Block diagram representation of the greenhouse control system.The system is composed of controllerS M , energy management system S E , greenhouse climate and lighting system model S G and crop growth and transpiration model S C .The control inputs to the greenhouse system are denoted by u.The elements of u input to S C encompass the harvest of fruits and leaves.The inputs to the greenhouse climate model, e.g.ventilation, screen deployment, heating, CO 2 injection are denoted by u g , these are internal to the model F. The uncontrollable inputs to the greenhouse climate model are denoted by d, the outside weather.Variables y g and y c denote the effect of the greenhouse climate on the crop (temperature, CO 2 concentration, radiation and relative humidity) and the effect of the crop on the greenhouse (assimilation and transpiration), respectively.Kuijpers et al. (2021).The increased update rate is a first step towards mitigation of the effect of the forecast error (Mayne, 2014).At every sample, the optimisation algorithm obtained feedback of the system through the state of the system at the present time x t .The state of the system includes the effect of the weather forecast error at previous time instances.The realised forecast error e2R n d , the value of which is uncertain, manifests itself through b dðjjkÞ ¼ dðj þkÞ þ eðjjkÞ.One of the approaches to handle the uncertainty is to neglect its effect and design a controller based on the nominal case, assuming dðj þkÞ ¼ b dðjjkÞ.The latter technique, often combined with a high update rate as discussed previously, is referred to as certainty equivalence MPC (CEMPC) (Saltık et al., 2018) and mitigates the effect of the forecast error eðjjkÞ on the operational return J through feedback.The CEMPC controller was used in the simulation studies as the approach which neglects the effect of the uncertainty on the performance of the controlled greenhouse system.
Throughout the results presented in this paper, three indicators are used to compare the performance of the greenhouse system under varying configurations and conditions.

The first indicator is the operational return
which balances the income through fruit yield and electricity sales with the cost for buying electricity, using gas and buying (pure) CO 2 .To isolate the effect of uncertainty in the weather forecasts, the uncertainty in price forecasts (e.g.c frt , c eby ) is omitted (Kuijpers, Antunes, et al., 2021).The latter indicator is also employed as objective function in the optimisation algorithm (2).The following two indicators will only be used to compare simulations of the greenhouse system and have been chosen in accordance with Kuijpers et al. (2021) to allow for comparison.The second indicator is the gas use denoted by sð,Þ : R nu /R (m3 m À2 ), gas is used by the CHP u chp and boiler where a g ¼ 31:65 MJ m À3 (Vermeulen, 2016) represents the energy content per cube of gas.The gas use is defined here as a separate performance indicator, but is also required for both lð ,Þ and p 1 ð ,Þ and is included in c chp and c boi .The third indicator is the carbon footprint p 2 ð ,Þ : R nu /R (kg:m À2 ) The carbon footprint increases with the use of gas, buying pure CO 2 u cby and buying electricity u eby , the carbon footprint decreases with selling electricity u ese .All aforementioned indicators are linear combinations of the inputs.The weights c used to express the contribution of each input to indicator value in the respective unit are defined in Kuijpers et al. (2021) and Vermeulen (2016, p. 330).

Weather forecasts & observations
In this subsection, the origin of the historical data which were used in the estimation of parameters of the stochastic weather forecast error model is presented.The historical data for the weather forecasts originate from the HARMONIE-AROME Cy40 model 3 run by the KNMI.The weather forecasts are published every 6 h starting at midnight, i.e. at time instances k2f0; 6M; 12M; …g, where M ¼ 4 is the number of time instances in an hour.The lead time of a forecast t f 2R represents the time interval between the time instance the The grid consists out of 340 Â 340 points spaced 0:05 along both lateral and longitudinal axis.Although the research presented here involves a simulation study and the greenhouse is not location-bound, we chose the grid point closest to the greenhouse where the data from the weather realisation d originates from, Bleiswijk, grid point 52 2 0 9:6 00 N; 4 30 0 50:36 00 E (Kempkes et al., 2014).The historical data for the weather observations originated from the automatic weather station Rotterdam (06344)4 located at 51 58 0 N; 04 27 0 E, 9 km away from the forecast grid point.All observations and 400 forecasts, here denoted by d2R n d , for that specific location were stored in the period from 01 À 01 À 2020 to 01 À 07 À 2020.In Fig. 2, a subset of this data, i.e. historical temperature observations and 33 temperature forecasts, is presented.Every individual forecast is represented by a line of 48 h in length.The historical data for the forecasts and the observations were used to identify the stochastic model of the weather forecast error as described in Subsection 2.3.
In the simulation studies, the realisation of the weather applied to the system d originated from an experiment described in Kempkes et al. (2014), where various energysaving options in greenhouses were investigated in a Venlow Energy kas located in Bleiswijk, The Netherlands.The data, which consist of outside air temperature T out ð CÞ, outside absolute air humidity H out ðg:m À3 Þ, outside air CO 2 concentration CO 2;out ðg:m À3 Þ, wind speed v wind ðm:s À1 Þ and global radiation Q sun (W m À2 ), are measured at 5 min interval, during the years 2011 to 2014.In the simulation studies here, only year 2014 was used.The historical observations by the KNMI could not be used for this purpose due to the unavailability of CO 2 concentration measurements and an insufficient temporal resolution, i.e. the observations by the KNMI have a resolution of 1 h where at least 15 min is required.Subsection 2.3 details how the data from Kempkes et al. (2014) was combined with the stochastic weather forecast error model to synthesise weather forecasts for the realisation of the weather.
In this study the forecast data from the HARMONIE-AROME Cy40 model between 01 À 01 À 2020 and 01 À 07 À 2020 was used.This not only limits the data to the meteorological winter, spring and summer seasons but also does not capture the year-to-year differences in the climate or the forecasting of it.

Weather forecast error model
For a weather forecast published at time instance k and forecasting weather variables at time instance k þ t f , the forecast error e is defined as the difference between the published weather forecast d and the corresponding weather realisation d, i.e.
The forecast error for t f > 0 can only be calculated in retrospect as dðk þt f Þ is not available at time instance k.A stochastic forecast error model was used to describe the stochastic properties of the error in the weather forecast.The structure of the stochastic forecast error model was chosen to reflect that of an autoregressive model of order 1 (ARð1Þ), a model with a similar structure was successfully employed for building climate control in Oldewurtel et al. (2014) and Zhang et al. (2013).The deterministic part of this ARð1Þ-model is able to capture structural deviations such as a non-collocated forecasting and measuring facility (Doeswijk, 2007).The latter three studies have shown the ability of an ARð1Þ-model to describe the weather forecast error.The ARð1Þ-model is where wðt f Þ $ N ðm h ; s h ðt f ÞÞ represents the residual after parameter estimation with mean m h 2R n d and standard deviation s h 2R n d .The contribution of this research is the estimation of the parameters J, x and wðt f Þ using relevant historical observations and forecasts from a Dutch weather forecasting service.With respect to the model proposed in Oldewurtel et al. (2014) and Zhang et al. (2013) the dependency of w on t f was added to account for the increase in the forecast error e for increasing values of the lead time t f .As this paper describes a simulation study that aims to quantify the effects of forecast errors on a controlled greenhouse system, an ARð1Þ-model will suffice in removing unwanted, unrealistic effects such as from a non-collocated forecasting and measuring facility.The estimation of the parameters J, x and wðt f Þ in ( 7) can be cast into a least squares regression problem and is solved as such.
To support the hypothesis that an ARð1Þ model structure suffices to describe the correlation between the data samples, the partial autocorrelation function 5 (PACF) was calculated for the available historical data.The model in (7) describes the evolution of the forecast error for T out , absolute air humidity H out , wind speed v wind and global radiation Q sun , all elements in d except for the CO 2 concentration CO 2;out as no forecasts of the CO 2 concentration are available in the historical forecast data.The left panel of Fig. 3 depicts the PACF for the temperature element in e, the lags follow from the lead time t f in terms of hours, as the temporal resolution of the forecasts is 1 h.The PACF for all elements in e is significant for lag 1 only, hence an ARð1Þ model structure suffices to describe the correlation between the samples in e (Box et al., 2016).The PACF of the temperature component of w in the right panel of Fig. 3 shows that partial autocorrelation at lag 1 is compensated by the ARð1Þ-model.This supports the hypothesis that the parameters were estimated successfully and that the model structure is sufficient to describe the autocorrelation in the data.The PACF of the residuals also supports the choice of model order, first, and the fact that no moving average mechanism is present.Results leading to similar conclusions were obtained for the other components of vectors e and w.The interactions between weather variables are not captured by the individual PACFs, therefore J in (7) was chosen to be diagonal.Also, the residuals w were modelled as independent signals.An analysis, including the cross terms in the correlation of e, was out of scope of this research as the model in (7) provides a description of the weather forecast error which is sufficiently accurate for this research.After identification, the standard deviation s h of w is determined per value of t f .Propagation of the stochastic model allows for the creation of artificial forecasts d with similar stochastic properties as the historical weather forecasts based on a weather realisation d, i.e. d ¼ d þ e and substituting e with (7).In order to generate a forecast, realisations of wðt f Þ $ N ðm h ; s h ðt f ÞÞ have to be obtained.In Fig. 4, the historical temperature observations and 33 synthesised forecasts are presented for the same period of time as in Fig. 2. Every individual forecast is represented by a line of 48 h in length.
In order to compare the stochastic properties of the historical forecasts and the synthesised forecasts, Fig. 5 presents the mean m and corresponding 3s-bounds as a function of t f for both forecasts.The left panel of Fig. 5 presents the average forecast error, denoted by m h , and 3s-bounds of for the temperature element in the historical forecasts.The right panel of Fig. 5 presents the same stochastic properties of the synthesised temperature forecast errors.One can observe that the stochastic properties of the synthesised Fig. 3 e Partial autocorrelation function (PACF) of (left) historical forecast error e, (right) the residual after parameter estimation w.The lags follow from the lead time t f and are in terms of hours.From the left panel one can observe that an ARð1Þ-model can sufficiently capture the correlations between different lags in eð,jkÞ.The right panel shows that the residual after parameter estimation w is not autocorrelated.The PACF presented here is for the temperature element, the PACF for the other weather variables have not been plotted. 5The partial autocorrelation function at lag k2N of variable zðtÞ2R with t2N is defined as the correlation between zðtÞ and zðt ÀkÞ with the correlation between zðtÞ and zðt À1Þ; zðt À2Þ; …; zðt Àk þ1Þ removed (Box et al., 2016).temperature forecast errors are similar to those of the original forecasts.

Kalman filter
The parameter estimation in Section 2.3 shows that a part of the weather forecast error exhibits deterministic behaviour, e.g. because of non-collocated measurements or deterministic errors within the weather forecast itself.The Kalman filter presented in this section is able to compensate for this deterministic part.The Kalman filter provides a better estimate than one could obtain without, i.e. by simply relying on the measurements.The previously published forecast synchronised with the prediction horizon is denoted by dðjjkÞ2R n d and defined as dðjjkÞ where a2N is the number of time instances since the previously published forecast.P2R np,n d Ânp,n d and P2 R np,n d Ânp,n d , respectively.The Kalman filter is graphically represented as part of the block diagram in Fig. 6.The switch in Fig. 6 represents the ability to exclude the effect of the Kalman filter in the forecast.The left part of Fig. 6 represents the synthesis of artificial forecasts as outlined in Subsection 2.3.
To arrive at the error model for the prediction stage of the Kalman filter, the model presented in (7) was augmented with all values for t f .The Kalman filter iterates along the time axis k instead of the lead time axis t f .The forecast error at time instance k, eðt f kÞ, however, can be related to eðt f k þ1Þ by substitution of eðt f À1 k þ1Þ ¼ eðt f kÞ in (7).These operations, concerning both the t f -timescale and the k-timescale, are depicted in Fig. 8. Figure 8 depicts the transition from eðΤjKÞ (indicated by 1) to eðΤjK þ1Þ (indicated by 2), using equation ( 7) and eðt f À 1 k þ1Þ ¼ eðt f tÞ.After the two operations, the prediction model used in the prediction stage is obtained By assuming that the realised forecast error eðt f tÞ results from the same ARð1Þ process as in ( 7), the evolution of the estimation error of the forecast error, eðt f kÞ À b eðt f kÞ; can be described as which is an ARð1Þ-model.The covariance of a forecast by an ARð1Þ-model is given by b according to Box et al. (2016).The elements in the estimate covariance matrix b Pð0j0Þ in ( 12) are given in ( 14), in which j and k follow from the multiplication of two vectors having elements similar to the vector presented in (9).The process covariance matrix of the Kalman filter represents the covariance of w, i.e.Q ¼ ½s h ð0Þ; s h ð1Þ;…;s h ðn p Þ.The elements of w are zero-mean Gaussian processes as required to obtain optimal state estimation (Anderson & Moore, 2005).The measurement covariance matrix of the Kalman filter was chosen as R ¼ 0:1 5 I np , representing the accuracy of the climate sensors.In the simulation studies, however, the measurements were equal to the actual realisation of the climate, as indicated in Fig. 7.

Randomised MPC
The RMPC controller was used in the simulation studies as the approach which takes the effect of the uncertainty on the performance of the controlled greenhouse system into account.The RMPC algorithm is presented in Zhang et al. (2013) and discussed in Saltık et al. (2018) and Schildbach et al. (2012).This subsection details the application of the RMPC approach to the greenhouse control problem, a more elaborate description of the RMPC algorithm is given in Zhang et al. (2013).The RMPC algorithm handles the uncertainty in the optimisation problem by drawing N sc independent and identically distributed (i.i.d.) samples from the full-horizon uncertainty space, i.e. number of scenarios, aims to optimise a (single set of) optimised trajectories for the controllable inputs of the system u * d .As the operational return does not depend on the state trajectories, the objective value for each scenario will be equal.The state trajectories x r d cr2f1; …; N sc g in each of the scenarios will be different, and the state trajectories of all scenarios should satisfy the constraints.The optimisation problem in (2) is rewritten to include multiple scenarios subject to: x r d ð0jkÞ ¼ x t cj ¼ f0; … ; Ng; cr2f1; …; N sc g: In literature, various robustness properties have been established for RMPC with respect to chance constraints.For example, Zhang et al. (2014) establish a relation between the number of scenarios N sc and the confidence in only violating a pre-specified number of constraints.An increased number of scenarios would increase the confidence that the constraints of the greenhouse system would not be exceeded.However, due to the scale of the optimisation problem in (15), induced by the prediction horizon, the model size and number of constraints, obtaining N sc > 3 is computationally not viable in our implementation.Therefore, instead of replacing our set of constraints hð ,Þ and the lower-and upper bounds of the design variables by chance constraints, the N sc ¼ 3 scenarios were required to satisfy all the constraints.Even though the full-horizon uncertainty space is not fully sampled by using N sc ¼ 3 scenarios, we hypothesise that the resulting optimised control trajectories from (15) are more robust to various forecast errors from the uncertainty space as compared to the trajectories resulting from (2).

Forecast types & simulation studies
In order to evaluate the effect of the weather forecast error on the performance of the controlled greenhouse system, various  configurations of the system were simulated.In this subsection, the various types of forecasts and simulation studies are presented.Four forecast types were defined that distinguish the various configurations of the system.The forecast is linked to the configuration of the controller, e.g. the use of multiple scenarios is linked to the RMPC algorithm, the combination of the two is, however, referred to as a forecast type in the remainder of this paper.The different forecast types are performance bound (PB), Certainty Equivalence MPC (CEMPC), CEMPC in combination with a Kalman filter (CEMPCþ KF) and RMPC in combination with Kalman filters (RMPC þ KF).The main features of the forecast types are presented in the first four columns of Table 3.Note that all forecasts are generated using the (7), the realised CO 2;out was substituted for its forecast.The forecast types provide the optimisation problem in (2) or ( 15) with different values for b dðjjkÞ: Performance Bound (PBMPC): b dðjjkÞ ¼ dðk þjÞ, the forecast of the weather b d matches the realisation d, resulting in a non-causal, and therefore not realistic, control approach.This forecast type, with full prior knowledge of the weather, resembles the maximal attainable performance.
Certainty Equivalence MPC (CEMPC): b dðjjkÞ ¼ dðjjkÞ, the previously published forecast synchronised with prediction horizon of the optimisation algorithm is assumed to be equal to the realisation d.The published forecasts are synthesised using (7).This forecast type resembles the case in which the effect of the uncertainty on the performance of the controlled system is neglected by assuming the forecast matches the realisation.The weather forecast error is a result of both the deterministic and stochastic part of the model in (7).CEMPC in combination with a Kalman filter (CEMPCþ KF): the previously published forecast synchronised with the prediction horizon is updated using the estimated forecast error by a Kalman filter, i.e.This forecast type takes into account the effect of the uncertainty on the performance of the system in the synthesis of the controller.
The Kalman filter is used to update the forecasts published every 6 h.Linear interpolation is used determine the quarterhourly values of the forecast error from the hourly values from the (updated) forecasts.The length of the KNMI forecasts was 48 h, the values for wðt f Þ therefore are estimated up to t f ¼ 48M.In using (4) to create forecasts with a length similar to the prediction horizon of (2), N ¼ 288 (72 h), the standard deviation of the stochastic part, i.e. s h was kept constant for j > 48M, i.e. s h ðt f ,MÞ ¼ s h ð48 ,MÞct f !48.The latter was schematically represented in Fig. 9.
Three simulation studies were used to obtain a better understanding of how the forecast error affects the controlled greenhouse system: simulation study 1: the stochastic properties of the original forecasts and the forecasts filtered by the Kalman filter are compared.simulation study 2: this study provides insight into the effect of the forecast error on the value optimised control trajectory u * d at the first time instance, i.e. u * d ð0jkÞ.The sensitivity of u * d ð0jkÞ to changes in the forecast error indicates to what extent the optimised control trajectories are affected by the forecast error.This simulation study is performed for both original forecasts CEMPC and forecasts filtered by the Kalman filter CEMPC þ KF. simulation study 3: the controlled greenhouse system was simulated during three 7-day intervals, for each of the four forecast types.The performance bound is used as reference to compare the performance of the controlled system to the other forecast types.In this simulation study, the performance bound represents the maximal attainable performance, resulting from a static optimisation based on the optimisation problem Table 3 e The various forecast types employed in this research, performance bound (PB), certainty equivalence MPC (CEMPC), CEMPC in combination with a Kalman filter (CEMPC þ KF) and randomised MPC in combination with Kalman filters (RMPC þ KF).The columns specify, the origin of the forecast supplied to the controller b d, whether it is updated, how it is interpolated from hourly to quarter-hourly values and how many scenarios are employed.The three last columns denote in which simulation study the forecast type is used (-)

Results
The results of the simulation studies introduced in Subsection 2.6 are presented in this section.The subsections are ordered accordingly.

Forecast error and estimation
The stochastic properties of the weather forecast error presented in Subsection 2.3, more specifically in Fig. 5, are valid at the time instances at which a forecast is published.In between those time instances, the previously published forecast was synchronised with the prediction horizon of the optimisation algorithm, according to (8).Also, if the Kalman filter was enabled, the previously published forecast was updated using local weather measurements.In the simulation study presented in this subsection, the stochastic properties of the forecast error are evaluated while taking into account all time instances and thus including the time instances at which no forecast was published.Figure 10 presents the stochastic properties, mean m and 3s-bounds, as a function of t f for (left) the unfiltered temperature forecasts (right) filtered temperature forecasts for all time instances.To evaluate the stochastic properties, the 2016 forecasts from the three 7-day simulations presented in Subsection 3.3, were evaluated.Through using the forecasts from the 7-day simulations, a realistic distribution was obtained between time instances at which a forecast is published and time instances at which no forecast is published.The effect of using the previously published forecast until a new forecast is published can be observed when comparing the left panel of Fig. 10 to the right panel in Fig. 5.One can observe from Fig. 10 that the uncertainty, in terms of 6s, exceeds 8 C at a lead time t f ¼ 31 (7.75 h).When evaluating only the time instances at which a forecast is published, see Fig. 5, the uncertainty is lower, i.e. 6s exceeds 8 C at a lead time t f ¼ 76 (19 h).The 3s bounds, thus, grow faster for increasing lead times t f when evaluating all time instances.The latter is expected as the lead time of the published forecast grows as a in (8) increases.The 3s-bounds increase and decrease periodically, an effect which is hypothesised to be due to the updates arriving every 6 h.The mean forecast error m is expected to reach a constant value of À 0:38 C, which can be found through equating eðt f þ1 kÞ to eðt f kÞ and solving for eðt f kÞ in (7).The constant value of the forecast error in the left panel of Fig. 10 for t f > 200, however, is À 1:09 C, significantly lower.How this links with the synchronization of the published forecast with the prediction horizon of the controller has not been the subject of further study.
The Kalman filter aims to improve the forecasts through (a) updating the previously published forecasts using local weather measurements and (b) compensating the effect of the deterministic part in (7).In Fig. 10, one can observe that the contribution of the deterministic part in ( 7) is relatively small compared to the uncertainty, i.e. the mean deviates in the order of 1 C whereas the contribution of the stochastic part, measured in terms of 3s, is in the order of 4 C.The second aim of the Kalman filter is, therefore, expected to contribute less to the overall improvement.One can observe the effect of the improvements (a) and (b) when comparing the left panel of Fig. 10 to the right panel.The updates of the forecasts decrease s for t f ¼ 0 from 0:96 C to 0:49 C. The effect of compensating the deterministic part of the stochastic weather forecast error model can be observed by comparing the average forecast error, indicated by m, its value for t f > 200 decreased from À1:09 C to À 0:71 C.

Uncertainty analysis
The receding horizon optimal controller based on ( 2) is sensitive to the forecast error if the value of u * d at the first time instance is considerably different for various realisations of the forecast error.This sensitivity was evaluated by repeatedly solving the optimisation algorithm in (2) for different values of the forecast error.As only the value of u * d at the first time instance, i.e. u * d ð0jkÞ is applied to the system, this analysis was limited to u * d ð0jkÞ.The sensitivity with respect to the operational return is not included as this cannot be interpreted in a simple way, because the operational return is the result of a closed-loop system.Moreover, analysing the sensitivity per control input can support reasoning about the sensitivity of specific subsystems within the total system.As the sensitivity may depend on the state of the system x t and the prevailing weather d (not the forecast error), the sensitivity is evaluated at distinct time instances throughout the year.Five time instances were selected, 17th of February at 00:00 h, 12th of March at 00:00 h and 12:00 h and the 12th of June at 00:00 h and 12:00 h.The state of the system x t at the various time instances was based on the HPS simulation presented in Kuijpers et al. (2021).For each of the time instances, 80 simulations were performed using forecasts updated by the Kalman filter, denoted by CEMPC þ KF, and the original forecasts denoted by CEMPC.Generally, forecast errors are thought to deteriorate the performance of the Fig. 9 e Schematic representation of the weather forecast used in the prediction of the RHOC.The (filtered) weather forecast is defined at hourly intervals, represented by the green instances above.In-between the estimates the hourly values are linearly interpolated to arrive at values for every 15 min.After 2 days (¼48M) instances into the prediction horizon, s h is kept constant at its value s h ð48 ,MÞ.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)b i o s y s t e m s e n g i n e e r i n g 2 1 4 ( 2 0 2 2 ) 2 0 7 e2 2 9 automatic climate control, we therefore expect differences in the control input.
From the elements in the vector u * d ð0jtÞ, see Table 2, only those relevant to understand the effect of the forecast errors are presented.Figure 11 presents histograms of the value of the level of operation of the CHP u chp , level of operation of the boiler u boi , electricity exchange with the grid u eby À u ese , CO 2 injection u CO2 , energy flux to the heat buffer u sto , ventilation rate u ven , screen deployment u scr and leaf harvest u lea in u * d ð0jtÞ.If a specific element of u * d ð0jtÞ attains a comparable value for the various forecasts, i.e. the element is not sensitive to forecast errors, a narrow distribution results.If a specific element of u * d ð0jtÞ attains various values for various forecasts, i.e. the element is sensitive to forecast errors, various columns are shown with a combined sum of 80. Figure 11 presents histograms for simulations representing the state of the system and the prevailing weather at February 17, 00:00 h, similarly Figs. 12 and 13 present histograms for March 12, 12:00 h and June 12, 00:00 h, respectively.The width of the histogram bins result from the minimal and maximal bounds in Table 2, except for the ventilation rate u ven as the minimal and maximal ventilation rate depend on the windspeed.Figures 11e13 show the differences in u * d ð0jtÞ for forecasts updated by the Kalman filter, denoted by CEMPC þ KF, and the original forecasts denoted by CEMPC.The histograms of the simulations with the filtered forecasts are not considerably different from the histograms of the simulations with unfiltered forecasts.One can thus conclude that the difference in sensitivity to forecast error between the two forecast types is negligible.The subsequent results are therefore independent of the forecast type.
One can observe from Figs. 11e13 that, overall, the elements such as the level of operation of the CHP u chp , boiler u boi and CO 2 injection u CO2 were sensitive to the forecast error.We hypothesise that these are sensitive as these affect the greenhouse air temperature.The sensitivity of the CO 2 injection u CO2 is due to the dependency of CO 2 concentration on the greenhouse air temperature through the ideal gas law.Note that no forecast error is present in the forecast of the outside air CO 2 concentration.We hypothesise that the sensitivity of CO 2 injection u CO2 with respect to the forecast error is due to the CO 2 concentration being operated close to its maximum value, see Table 1, a different temperature therefore requires different u CO2 to meet the constraints.Although during most time instances, all values of the inputs elements were in close proximity of the simulation with the weather realisation d, considerably different values were obtained for especially for the level of operation of the CHP u chp , and CO 2 injection u CO2 .
Figures 15 and 16 in Appendix A show, similarly to Fig. 11, the uncertainty analysis for March 12, 00:00 h and June 12, 12:00 h.Overall, one can observe an increased sensitivity of the elements in u * d ð0jtÞ to the forecast error as the season progresses.We hypothesise that this is due to the system being operated closer to constraints, as discussed previously.Throughout the season, the system is operated closer to the relative humidity bound, the data on this, however, is ambiguous as this bound is not active in all cases in which considerably different trajectories were observed.

Controlled system performance
The simulation study presented in this subsection evaluated the effect of weather forecast error on the performance of the controlled greenhouse system.As the effect of forecast errors on the performance of the system depends on the state of the system x t and the prevailing weather d, three intervals of 7 days spread throughout the growing season are presented.The three intervals are the 11th up to and including the 17th day of the months January, April and May.The state of the greenhouse and the crop at the start of the interval were based on the simulations with HPS lighting in Kuijpers et al. (2021).To prevent effects of the forecast error after the 7 À day interval, the horizon shrank to ensure no time instances after the 7 À day interval were considered in the prediction horizion.
A successful execution of the optimisation algorithm in (2) guarantees constraint satisfaction of the system when the resulting trajectory is applied and the prevailing weather is Fig. 10 e The stochastic properties, mean m and 3s-bounds, as a function of t f for (left) the unfiltered temperature forecasts (right) filtered temperature forecasts, for all time instances.The stochastic properties in these panels are based on 2688 forecasts each.When compared with the left panel, the right panel shows the effect of the Kalman filter.
equal to its forecast.The control approaches without perfect weather predictions, however, may potentially cause constraint violation.Therefore, after optimisation, the first element of the optimised control trajectory, alongside the realisation of the uncontrollable input, were used in an openloop simulation to determine the state of the system at the next time instant.Table 5 presents the integrated constraint violation values for the three weather forecast types without perfect predictions, i.e.CEMPC, CEMPC þ KF and RMPC þ KF, for the three 7 day intervals.One can observe that the upper bound of the state constraint on CO 2;air is violated during all of the three 7 day intervals, the violation values itself, e.g.5:04 g:h: m À3 for CEMPC forecast type in April are small compared to the magnitude of the variable CO 2;air .The violations in the lower Fig. 11 e Histograms of the value of various elements in u * d ð0jtÞ, for various realisations of the forecast error.The prevailing weather and state of the system used in these simulations represent the actual system at February 17, 00:00 h.For both the unfiltered forecasts CEMPC and the filtered forecasts CEMPC þ KF, 80 simulations have been performed.The value of the (left-to-right, top-to-bottom) level of operation of the CHP u chp , level of operation of the boiler u boi , electricity exchange with the grid u eby À u ese , CO 2 injection u CO2 , energy flux to the heat buffer u sto , ventilation rate u ven , screen deployment u scr and leaf harvest u lea are presented.The value at the first time instance of the respective elements for the simulation with zero forecast error is represented by the orange line.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)Fig. 12 e Histograms of the value of various elements in u * d ð0jtÞ, for various realisations of the forecast error.The prevailing weather and state of the system used in these simulations represent the actual system at March 12, 12:00 h.For both the unfiltered forecasts CEMPC and the filtered forecasts CEMPC þ KF, 80 simulations have been performed.The value of the (left-to-right, top-to-bottom) level of operation of the CHP u chp , level of operation of the boiler u boi , electricity exchange with the grid u eby À u ese , CO 2 injection u CO2 , energy flux to the heat buffer u sto , ventilation rate u ven , screen deployment u scr and leaf harvest u lea are presented.The value at the first time instance of the respective elements for the simulation with zero forecast error is represented by the orange line.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)bound of C frt represent the effect of impossible harvesting, i.e. harvesting when the fruit buffer C frt already depleted.As C frt is expressed in dry-mass, a violation of 0:002 kg m À2 represents 0:035 kg m À2 in fresh weight, which, in the case of January, results in 0:04 V m À2 income through yield.Therefore, a violation of the state constraints for C frt significantly affects the operational return, which, for a 7 day-interval in January is 1:16 V m À2 .Because of the magnitude of this effect on the operational return and the fact that the contribution of fruit harvest to the operational return is straightforward, the resulting operational return was compensated for the constraint violation.A violation of the state constraint on CO 2;air also affects the objective function, its effect on the objective function, however, is hypothesised to be smaller and more difficult to evaluate than the fruit harvest violation.
The forecast error combined with the prevailing weather may result in a non-realistic weather forecast which makes the optimisation problem in (2) infeasible.If the solver of (2) did not converge to a feasible solution, the optimised control trajectory u * d from the previous time step was used for the corresponding time step, e.g.u * d ðt þ1Þ for the first non-converged solution.The largest number of subsequent time steps at which the solver did not converge was 11, which represents an interval of 2:75 h of the aforementioned open-loop control strategy.The highest number of time steps within an interval, not necessarily subsequent, was 46 for the CEMPC forecast type in May, which is 7% of the time steps within the interval.
Table 4 presents the resulting operational return J as in (1), gas use s, carbon footprint p 2 , the integrated value of the ventilation rate over the interval and averaged screen use for the different forecast types, PBMPC, CEMPC, CEMPCþ KF, RMPC þ KF and time periods throughout the growing season.The operational return J in Table 4 is compensated for the violations of the lower-bound of C frt , based on Table 5.The simulations with the controllers without explicitly taking into account the uncertainty in the weather forecast error, i.e.CEMPC and CEMPC þ KF result in a loss of performance of at most 0:03 V m À2 with respect to the performance bound.Concluding, the effect of weather forecast errors on the performance of the greenhouse is small.The increase in the operational return of the simulations with CEMPC and CEMPC þ KF controller with respect to the performance bound in April can be explained through constraint violation, see Table 5.In practice, the operational return for CEMPC and CEMPC þ KF will be lower.Figure 14 presents optimised trajectories during the interval of January 11th to January 17th for the four forecast types simulated: PBMPC, CEMPC, CEMPC þ KF and RMPC þ KF.One can observe the trajectories coincide to a large extent, the most notable differences are during the night.We hypothesise that this is due to the low outside air temperature in the period, this will result in operation close to the lower temperature bound and which can be violated upon selection of different type of forecast.The spikes that can be observed in (among others) the fruit harvest u frt occur at times when a new forecast arrives from the KNMI.At these times the new forecast differed considerable from the forecast that was updated, leading to a different control strategy.
A detailed analysis of the data presented in Table 4 shows that the difference between original forecasts and updated forecasts, between CEMPC and CEMPC þ KF, is small, yielding not significantly different values for all indicators except for the integrated ventilation sum u ven and gas use.Overall, the constraint violation values for the simulations with CEMPC þ KF are lower than those for CEMPC, see Table 5.The latter indicates that the weather forecasts in the CEMPC þ KF are Fig. 13 e Histograms of the value of various elements in u * d ð0jtÞ, for various realisations of the forecast error.The prevailing weather and state of the system used in these simulations represent the actual system at June 12, 00:00 h.For both the unfiltered forecasts CEMPC and the filtered forecasts CEMPC þ KF, 80 simulations have been performed.The value of the (left-to-right, top-to-bottom) level of operation of the CHP u chp , level of operation of the boiler u boi , electricity exchange with the grid u eby À u ese , CO 2 injection u CO2 , energy flux to the heat buffer u sto , ventilation rate u ven , screen deployment u scr and leaf harvest u lea are presented.The value at the first time instance of the respective elements for the simulation with zero forecast error is represented by the orange line.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)closer to the actual realisation and thus result in a lower constraint violation value.
From Table 4, one can observe that the simulations with the RMPC controller result in a loss of performance of at most 0:04 V m À2 with respect to the performance bound.From Table 5, one can observe that the RMPC controller violates the least amount of constraints.The low amount of constraint violation is linked to the increased gas use and increased ventilation sum of the simulations with RMPC controller as the controller is using more resources to satisfy the constraints on the N sc scenarios.Concluding, the RMPC controller results in more conservative behaviour which shows through an increased gas use and ventilation rate.A comparison between the simulations presented in Section 3.3 with the PB and RMPC configurations is justified because of the low constraint violation values.

Discussion
Here the results presented in Section 3 are discussed to answer the main question of this research.In Subsection 4.1 the ability of the stochastic weather forecast error model to synthesise forecasts with similar properties as those from the KNMI is discussed.In Subsection 4.2, the performance of the Kalman filter is discussed.Subsection 4.3 discusses the main question of this paper, i.e. to what extent does the forecast error affect the performance of the controlled system.Additionally, Subsection 4.3 discusses the performance of an algorithm that explicitly includes the effect of the uncertainty in the controller synthesis.Subsection 4.4 discusses the implementation in practice of the various parts of the proposed approach presented in this paper.
Fig. 14 e From top to bottom, left to right, the operational return J, fruit harvest u frt , gas use s, electricity exchange with net u eby À u ese , ventilation rate u ven , screen use u scr , HPS lighting u hps and CO 2 injection u CO2 in the period January 13th to January 17th for the four forecast types simulated: PBMPC, CEMPC, CEMPC þ KF and RMPC þ KF.
b i o s y s t e m s e n g i n e e r i n g 2 1 4 ( 2 0 2 2 ) 2 0 7 e2 2 9 4.1.
Weather forecast error model Subsection 2.3 presents a stochastic weather forecast error model, in this subsection the implications of this model on the simulation studies are discussed.From Fig. 5, one can conclude that the stochastic properties of the weather forecast error model in (7) match the stochastic properties derived from the historical data from the KNMI.There exists, however, a significant difference between the original and synthesised forecasts, as one can observe by comparing the original forecasts in Fig. 2 with synthesised forecasts in Fig. 4. The spread in forecasts for a specific time instance that have been synthesised is larger compared to the spread in the KNMI forecasts.This is the result of the inability of the model in ( 7) to correlate various forecasts for the same time instance.To illustrate this, consider the temperature at the 11th of February at 15.00 h, in the 2 days before that moment, 8 forecasts are published which predict the temperature at that specific instant.From Fig. 2, one can observe that the 8 KNMI forecasts at the specific time instant do not vary considerably, as opposed to the synthesised forecasts, see Fig. 4.
The weather forecast error model does not correlate the forecast error eðt f tÞ and eðt f Àa t þaÞ where a2f0; 6M; 12M; …g.This deficiency of the model in (4) does not affect the second simulation study presented in Subsection 3.2, as no subsequent forecasts are employed there.The third simulation study is affected by the latter deficiency as it involves simulations over an interval.We hypothesise that the loss in performance due to weather forecast errors presented was overestimated.It is, however, unclear to what extent the loss in performance is overestimated.The magnitude of the overestimation is hypothesised to be small as the sensitivity to forecast errors of terms that affect the operational return is low, as presented in Subsection 3.2.The loss in performance when using more realistic forecasts, such as the real forecasts, is thus hypothesised to be lower.

Kalman filtering
The Kalman filter presented in Subsection 2.4 aims to improve the forecasts through (a) updating the previously published forecasts using local weather measurements and (b) compensating the effect of the deterministic part in (7).In Fig. 10, one can observe that the standard deviation of the forecast error for t f ¼ 0 is decreased from 0:96 C to 0:49 C.
Table 4 e The (left to right) integrated values of ventilation rate over the interval u i ven , screen set u scr , operational return J, gas use s and carbon footprint p 2 for the different forecast types, PBMPC, CEMPC, CEMPCþ KF, RMPC þ KF, and time periods throughout the growing season.The values indicate a loss in performance of maximal 0:03 V m À2 with respect to the performance bound and an inconsiderable difference between the CEMPC and CEMPC þ KF simulations.
Table 5 e Constraint violation values for the three weather forecast types without perfect predictions, i.e.CEMPC, CEMPCþ KF and RMPC þ KF, for the three 7-day intervals.The values presented in this table represent the constraint violation values integrated over the 7-day interval.The columns, indicated either by ↑ or ↓, refer to constraint violations of the upper or lower bound, respectively.During the interval, of length 168 h, both the lower-and upper-bound can be violated..Only constrains which are violated during the simulations are presented here, i.e., the CO 2 concentration CO 2;air , the leaf buffer C leaf , the fruit buffer C frt and relative humidity RH.The decrease of the standard deviation for small values of t f was hypothesised to improve the performance of the controlled greenhouse system.The uncertainty analyses presented in Figs.11e13, however, show no considerable difference between the simulations with CEMPC and CEMPCþ KF.Also, the simulation results in Table 4 do not show a considerable difference between the simulations using CEMPC and CEMPC þ KF, for the cases that were tested.Although the trajectories optimised by the CEMPC þ KF approach seem to violate constraints to a lesser extent, the data over the various simulations is not unambiguous.The latter unambiguity is due to quality of the estimations by the Kalman filter, this quality varies over the simulations due to the stochastic part of the model in ( 7).
The Kalman filter also compensates the effect of the deterministic part in (7).The mean forecast error over various forecasts is EðeÞ ¼ x,ð1 À JÞ À1 , which can be derived from (7) using EðwÞ ¼ 0. The mean forecast error over various forecasts is represented by m in Fig. 10.By comparing both panels in Fig. 10, one can observe that the mean forecast error for t f < 150 is compensated, but that the forecast error for t f !150 is still considerable.The latter effect is hypothesised to be due to the measurement covariance matrix R, the value of which may be chosen even lower to reflect a higher accuracy of the local weather measurements.Also, the mean forecast error m deviates in the order of 1 C whereas the contribution of the stochastic part is in the order of 4 C.The influence of the mean forecast error was marginal compared to the stochastic element.
One of the reasons that the effect of the Kalman filter is low is because of the quality of the forecasts.The forecasts by the KNMI have a low mean forecast error and a low forecast error for low values of the lead time.This is, partly, due to the small distance between the point for which the historical forecasts were provided and the point at which the historical observations were obtained.The effect of the Kalman filter will be larger if the forecasts are provided for a location far from the greenhouse, see Doeswijk (2007).The distance between the measurement station collecting the historical observations and the grid point in the weather forecasts was approximately 9 km.An increased distance between the point of observation and the point for which the forecasts are provided will decrease the accuracy of the forecasts, and potentially increase the mean forecast error (Doeswijk, 2007).

Performance loss
The third simulation study, presented in Subsection 3.3, shows that the loss in performance, in terms of operational return, for a non-zero weather forecast error is small.For the simulations with CEMPC and CEMPC þ KF, the performance loss is at most 0:02 V:m À2 with respect to the performance bound, as can be observed from Table 4.The performance loss in practice is, however, hypothesised to be larger due to constraint violations, see Table 5.This, however, does not imply that the effect of the forecast error on the greenhouse system is small.The CEMPC controller is characterised by a high update rate (15 min À1 ).The effect of the uncertainty in the forecast error can, thus, be measured at the states of the system after 15 min, which allows the receding horizon algorithm to mitigate the effect of the uncertainty through feedback.The performance of a system controlled with a lower update rate will be lower than the performance specified here.Also, the feedback from the system, via x t in (2), does not include any error which is unrealistic.
The simulation studies presented here are based on forecasts from the KNMI, which are sufficiently accurate for this purpose, as mentioned the weather is forecasted at a location 9 km from the weather observation station.Also, the KNMI updates the weather forecasts every 6 h with the new output of the numerical weather model.
From Table 4, one can observe that for the simulations with RMPC þ KF a performance loss of at most 0:04 V m À2 with respect to the performance bound was obtained.By including multiple scenarios the controller becomes more conservative (Saltık et al., 2018) as the constraints in (15) have to be satisfied for all the N sc scenarios.The latter resulted in a higher ventilation sum and higher gas use throughout the interval for the RMPC þ KF controller as can be observed in Table 4.The conservatism induced by the scenarios considered in the RMPC þ KF controller did not significantly affect the operational return.Due to the increased gas use, the carbon footprint did increase significantly.The optimised control trajectories of the RMPC þ KF controller ensure a lower degree of constraint violation as can be observed in Table 5.Because of the diverse degrees of constraint violation for the simulations, the performance of the CEMPC and RMPC simulations should not be compared directly.
The uncertainty analyses show a considerable sensitivity for some of the elements in the first time instance of the optimised control trajectories, the resulting effect on the performance of the controlled greenhouse system is low.The latter observation is a result of the chosen objective function which is based on a weighting of the inputs according to l in (1).For example, in the case of Fig. 11, CO 2 injection u CO2 varies throughout the region constraint by U, i.e. highly sensitive.u CO2 affects J through crop growth and thus through u frt , this effect is small compared to e.g. the cost of gas.A high sensitivity to the forecast error might result in frequently changing control strategy if the forecast error is also varying (e.g. through updates by the weather service, updates from a Kalman filter).In the case of, e.g., the level of operation of the CHP u chp this might not be desirable from a maintenance perspective, as switching the CHP on and off might affect its durability (van Beveren et al., 2019).The latter is an additional reason to opt for a RMPC þ KF controller as its trajectories will be, generally, less sensitive to changes in one of the scenarios.
The model of the system that was used in this research describes the operation of the greenhouse with only fast time scales (Kuijpers et al., 2021).If slow time scales are also included, for example through the inclusion of a long-term energy storage or the development process of the crop, a longer prediction horizon will be required, e.g.weeks (van Straten & van Henten, 2010).A longer prediction horizon will require long-term weather forecasts, which typically are more uncertain as compared to the short-term weather forecasts used here.
The simulation studies all assume that the algorithms solving (2) and (15) find the global optimum for u * d ð0jtÞ, however due to the complexity of the problem this cannot be guaranteed.The settings of the NLP (and underlying LP) solver have been chosen with solely with the aim to find this global optimum (compared to e.g.speed of the algorithm).The global optimum and with that the results of the sensitivity analysis might shift when formulating a different control problem, e.g., choosing a different objective function.

Implementation
When implementing the approach presented in this paper, the weather forecasts will arrive from a weather forecasting service such as the KNMI.In this paper, however, a stochastic weather forecast error model was employed to (a) be able to generate forecasts based on a set of measured weather and (b) to extend this to multiple forecasts generated on one set of measured weather for the RMPC approach.
When employing an approach that neglects the uncertainty in the forecast error and only requires one forecast, i.e.CEMPC or CEMPC þ KF, one can use the weather forecasts arriving from a weather forecasting service.When the RMPC approach is employed, the scenarios might result from the various perturbations used in an ensemble prediction system such as presented in Frogner et al. (2019).Both approaches will remove the need for a stochastic weather forecast error model and the corresponding downsides mentioned in Subsection 4.1.
The RMPC approach in this research used N sc ¼ 3 scenarios, the inclusion of additional scenarios was computationally not viable in our implementation.With this configuration, a less constraint violations were observed, see Subsection 4.3.Further research on robust performance of the greenhouse system, e.g. through including more scenarios or including other control approaches (Saltık et al., 2018), should be performed to support these observations.
The receding horizon optimal control problems in (2) and (15) assume infinite computational power, as the state information x t arrives at the same time instant at which the first part of the resulting control trajectory u * d ð0jkÞ is applied to the system.In practice, the latter will induce a delay in the control system which will affect the performance of the system.A similar effect can be observed in the observed in the Kalman filter where the local weather measurement is used to update the weather forecast.
The results presented here apply to the greenhouse design and the climate type employed in this study.The parameters used for the greenhouse climate mainly originate from the model presented in van Beveren et al. (2015).The crop parameters originate from Vanthoor (2011) and limit the results of this study to tomato crops.The results presented are also based on the Dutch climate, which is classified by the K€ oppeneGeiger climate classification system as a temperate climate without a dry season and warm summers (Beck et al., 2018).The description of the forecast error and its stochastic properties is different for other locations with other climate types or when different weather forecasting service are used.The ideas proposed in this paper, however, may transfer to climates with a different classification as well.There has been no other study using a similar approach in different climates.As the models are white-box models, one could integrate different crops or greenhouses into this approach.

Conclusion & recommendations
The aim of this research was to quantify the loss in performance of the system due to weather forecast errors.Through simulation of three 7 day-intervals, spread throughout the growing season, we observed a loss in performance of at most 2 % with respect to maximal attainable performance in the optimally controlled greenhouse due to weather forecast errors.The 15 min update rate of the receding horizon optimal controller, in combination with the low forecast error in the forecasts by the Royal Netherlands Meteorological Institute contribute significantly to maintaining performance close to the maximal attainable performance.Due to the inability of the weather forecast error model to correlate various forecast for the same time instance, the employed forecasts are hypothesised to result in an overestimation of the loss in performance.The model of the greenhouse system employed in this research does not include slow time scales, hence a short prediction horizon is sufficient.This, however, considerably limits the required length of the weather forecast and reduces the magnitude of the uncertainty.The Kalman filter employed to improve the forecasts did not contribute considerably to the performance of the closed loop as compared to a configuration without Kalman filter.The Kalman filters contribution is hypothesised to be larger for weather forecasts with a lower quality, i.e., an increased mean forecast error.
The availability of weather forecasting service, with ensemble predictions, would remove the need for the stochastic weather forecast error model and corresponding assumptions.A next step would therefore be to apply RHOC in practice with real weather forecasts.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 e
Fig. 2 e Historical temperature observations d compared to 33 temperature forecasts d in the interval February 1st, 2020 to February 16th, 2020, originating from the HARMONIE-AROME Cy40 model.The forecasts have been plotted as separate lines starting at the time of the forecast with a length of 48 h.The first forecast d (light blue) is published at February 2nd at 06:00 h, the last forecast in this interval (pink) is published at February 14th at 00:00 h.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

TheFig. 5 e
Fig. 4 e Historical temperature observations d compared to 33 synthesised temperature forecasts d in the interval February 1st, 2020 to February 16th, 2020.The forecasts are plotted as separate lines starting at the time of the forecast with a length of 48 h.The first forecast d (light blue) is published at February 2nd at 06:00 h, the last forecast in this interval (pink) is published at February 14th at 00:00 h.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) written in a compact form as b eðkÞ ¼ Aẽðk À 1Þ þ Bx þ w (11) with augmented system matrix A ¼ I np 5J and augmented input matrix B ¼ I np Â1 5 I n d , where 5 denotes the Kronecker product and I np an identity matrix of size n p Â n p .In (11), w represents the augmented residual term of the ARð1Þ-model.The initial state estimate of the Kalman filter b eð0Þ was determined by propagating b eð0j0Þ ¼ dð0j0Þ À dð0Þ along t f using the ARð1Þ-model in (7).Let e2R np,n d represent the augmented realised forecast error.The initial value of the estimate Fig.6e Overview of various events that occur throughout time.The RHOC algorithm is run every t l ¼ 15 min.The forecasts specify the weather at hourly intervals Mt l ¼ 4t l , the Kalman filter operates at the same frequency and weather forecasts are published every 6 h.

Fig. 7 e
Fig. 7 e Block diagram of the synthesis of artificial forecasts d and the Kalman filter implementation.The block denoted by (4) uses the equation in (7) to generate artificial forecasts.The blocks denoted by [1] extract the first time instance from the time interval represented by the signal.The block denoted by (5) uses the equation in (8) to synchronise the previously published forecast d with the prediction horizon of the optimisation algorithm, resulting in d.The predicted and updated forecast errors are represented by b e and ẽ, respectively, their estimate covariance matrices by b P and P, respectively.The weather forecast input to the controller is denoted by b d.The weather forecast input to the controller b d can be chosen to be equal to the unfiltered forecast d, or equal to the updated forecast d À ẽ.

Fig. 8 e
Fig. 8 e An overview of the various operations on signals along the timescale k and timescale t f .The circles represent forecast instants, the arrows between them indicate the possible transitions.The AR(1)-model transforms forecasts instants along the t f -timescale (dashdotted line), the operation eðt f À1 k þ1Þ ¼ eðt f kÞ is represented by the dashed line.The Kalman filter operates along the timescale k. .
vmin ð ,Þ function describing the minimum ventilation rate as a function of wind speed, m 3 :m À2 :s À1 f vmax ð ,Þ function describing the maximum ventilation rate as a function of wind speed, m 3 :m À2 :s À1 Fð ,Þ discrete-time dynamical model of the greenhouse system R nx represents S E , S G and S

Table 1 e
States in the greenhouse system model and corresponding constraints represented by lower and upper bounds.States x d ¼ ½x s ; T air ; H air ; CO 2;air ; C buf ; C leaf ; C frt ; T c24

Table 2 e
Jong (1990)he greenhouse system model and corresponding constraints represented by lower and upper bounds.All bounds are fixed except for the bounds on the ventilation rate control input u ven .The bounds on u ven , f vmin ð ,Þ : R/ R and f vmax ð ,Þ : R/R, originate from the model by deJong (1990)and depend on the wind speed.Inputs u d ¼ ½u chp ; u boi ; u hps ; u led ; u eby ; u ese ; u CO2 ; u cby ; u sto ; u ven ; u scr ; u lea ; u frt forecast is published k 2 and the time instance at which the weather variable is predicted k 1 , i.e.t f ¼ k 1 À k 2 .The relevant weather variables forecasted are outside air temperature T out ð CÞ, outside absolute air humidity H out ðg:m À3 Þ, wind speed v wind ðm:s À1 Þ and global radiation Q sun (W m À2 ).The weather forecasts by the HARMONIE-AROME Cy40 model have a maximal lead time of 48 h, with a temporal resolution of 1 h, i.e. n p ¼ 49 time instances are forecast, at time instances t f 2 f0; 1M; 2M; …; 48Mg.The outside air CO 2 concentration CO 2;out ðg:m À3 Þ is not forecasted.The HARMONIE-AROME Cy40 model provides forecasts for a spatial grid covering Europe.
KF): N sc forecasts are generated using (7) and are updated by N sc independent Kalman filters.