Greenhouse air temperature predictive control using the particle swarm optimisation algorithm

The particle swarm optimisation algorithm is proposed as a new method to design a model-based predictive greenhouse air temperature controller subject to restrictions. Its performance is compared with the ones obtained by using genetic and sequential quadratic programming algorithms to solve the constrained optimisation air temperature control problem. Controller outputs are computed in order to optimise future behaviour of the greenhouse environment, regarding set-point tracking and minimisation of the control effort over a prediction horizon of 1 h with 1-min sampling period, for a greenhouse located in the north of Portugal. Since the controller must be able to predict the greenhouse environmental conditions over the specified time interval, it is necessary to use mathematical models that describe the greenhouse climate, as well as to predict the outside weather. These requirements are met by using auto regressive models with exogenous inputs and time series auto-regressive models to simulate the inside and outside climate conditions, respectively. These models have time variant parameters and so, recursive identification techniques are applied to estimate their values in real-time. The models employ data from the climate inside and outside the greenhouse, as well as from the control inputs. Simulations with the proposed methodology to design the model-based predictive air temperature controller are presented. The results indicate a better efficiency of the particle swarm optimisation algorithm as compared with the efficiencies obtained with a genetic algorithm and a sequential quadratic programming method. © 2005 Elsevier B.V. All rights reserved.


Introduction
Greenhouses are building structures that allow the creation of an indoor microclimate for crop development, protecting plants from adverse outdoor conditions. The greenhouse climate can be modified by artificial actuations, such as heating, ventilation and CO 2 supply, in order to provide the best environmental conditions. These modifications must be accomplished by spending additional energy in the production, requiring a regulator that minimises the energy consumption while keeping the state variables as close as possible to the optimum crop physiological reference. The use of model predictive control (MPC) for greenhouse indoor environment control has the advantage of providing the system with the ability to react before any deviations in the controlled variable take place, avoiding delays in the system response (Nielsen and Madsen, 1996). This class of control algorithms must employ models to describe and predict variables required for crop development over a specified time horizon. The MPC operation, within a process with bounded signals, usually involves the solution of a quadratic programming problem. This optimisation procedure is a fundamental part of model-based predictive control. The controller states are obtained by iterative numerical procedures that can be based on deterministic or stochastic algorithms. The optimiser must be able to handle constraints to model physical bounds, such as actuator saturation. Commonly, magnitude and rate constraints for the control actions and level constraints for the outputs are considered. Model predictive control cost functions, when subjected to restrictions, define a complex, non-linear, non-convex search space (Camacho and Bordons, 1994); hence, suitable for evolutionary algorithms optimisation. This paper reports the use of the particle swarm optimisation algorithm (Kennedy and Eberhart, 1995) to optimise a predictive controller within a greenhouse temperature control application. The results achieved with the proposed methodology are compared with the ones obtained with a genetic algorithm and a sequential quadratic programming method.

Experimental set-up
This section describes the experimental set-up and the methods used to simulate and control the air temperature inside a greenhouse located on the UTAD University Campus in the north of Portugal.
The greenhouse had a floor area of 210 m 2 and was covered with a 200 m polyethylene film. To control its climate, several actuators and sensors were installed and connected to an acquisition and control system based on a personal computer and a data acquisition card Advantech). This device had 16 analogue channels with 12-bit resolution, 31 digital input/output channels, one 16-bit counter, one analogue output and a real-time clock.
The weather and the greenhouse climate physical variables were acquired and/or controlled, with a sampling interval of 1 min. Inside and outside air and soil temperatures, relative humidities and solar radiation were measured using, respectively: PT100 sensors, capacitor sensors (F22H, Rotronic) and low-cost optoelectronic sensors (TSL230, Texas Instruments). Inside air CO 2 concentration was measured with an infrared gas analyser (GMP111, Vaisala) and the wind speed was measured with an anemometer (PV01-02, Electromatic).
Outside sensors were located 2 m above the soil. Inside the greenhouse, the air environmental data were measured at the central point 1 m above the soil, and soil temperature was measured at a depth of 5 cm. The data-acquisition system resolutions, including the sensors, were ±0.2 • C, ±2%, ±0.2 m/s and ±10 W/m 2 for temperatures, relative humidities, wind speed and solar radiation, respectively. The installed actuators were a ventilator with a flow rate of 38,000 m 3 /h, a shadow/thermal screen, a gas heating system with a heating power of 100,416 kJ/h, CO 2 injectors and an irrigation system.
Inside air temperature set-points, w(k), were constant during the daytime period and computed on-line for the night period. The desired air temperature for daytime conditions was 24 • C. For the night period, the desired air temperature depended on the outside air temperature, T o , and wind speed, W s , in order to reduce the heating system energy consumption. These values are computed each minute using Eq. (1), where W s (k) and T o (k) are the wind speed and outside air temperature measured at sample k.
The effective set-points used by the controller were obtained as the mean value of the last 10 calculations, which prevented high-frequency fluctuations in the air temperature caused by fast load variations in the wind speed or the transitions between night and day periods.

Air temperature model
To design model-based control algorithms to compute suitable control inputs in such a way that a cost performance function is optimised, it is necessary to use software models to compare different scenario outcomes. This requirement implies the use of dynamic models to predict the inside greenhouse climate, as well as the weather, in order to infer the optimal solution of a set of control sequences.
The dynamic changes in the greenhouse are determined by differences in energy and mass contents between inside and outside air, from exogenous energy as solar radiation or outdoor temperature and through the control actions taken. The greenhouse air energy balance is affected by energy supply and energy losses. The former is due to an artificial heating system and heat load imposed by the sun and the latter is due to transmission through the greenhouse cover and forced ventilation. Other energy and mass transport phenomena, for instance at the greenhouse soil, are neglected due to their small contribution to the overall air temperature dynamics.
Assuming that the greenhouse climate can be described by a linear system around an operating point, the greenhouse air temperature model is described by a parametric equation with exogenous inputs, u (Ljung, 1987), where T i is the air temperature, T the sampling interval, q −1 the backward-shift operator, n the number of delays from input to output, u the input signals vector and A and B are polynomials in q −1 . To select a model, a set of possible models, with different orders (degrees of polynomials A and B) and delays were first tested. There are several criteria that can be employed for this purpose, but the one most commonly used is Aikaike's Information Theoretic Criterion (AIC) (Akaike, 1974).
More explicitly, to simulate and predict the air temperature, T i , the model described by Eq. (3) is used, where a denotes the transfer function denominator parameters, b i the coefficients of the B i polynomials in the delay operator, is the difference between the heating pipe water (T p ) and inside air (T i ) temperatures, R o is the outside solar radiation, V and H the ventilation and heating control input signals, respectively, and the constant input 1 relates to the offset component. The model parameters a, b 1 , b 2 , b 3 and b 4 represent the partial contributions of each physical variable in the overall greenhouse air temperature and b 5 is the contribution of the continuous component. Since the model parameters are time varying (Boaventura Cunha et al., 2000), recursive identification techniques must be employed to estimate their values.
There are several methods that can be used to estimate the time-varying parameters of polynomials A and B (Åström and Wittenmark, 1995). For the case where the real parameters are slow time variant, as in this type of process, a recursive least-squares algorithm with exponential forgetting can be used (Boaventura Cunha et al., 1997). For the system described in Eq. (3), with delay n = 1, the model output,T i (k), can be expressed as: where: ϕ(k) is the data or regression vector and θ is the vector that contains the estimated parameters, The recursive computations of the vector parameters can be realized using Eqs. (7.1)-(7.3) known as Recursive Least-Squares algorithm (RLS), where K(k) denotes a gain matrix and P(k) is the covariance matrix of the estimated parameters. This RLS algorithm assumes that the process model parameters are constant. Since the real processes have time-varying parameters, the least-squares method could lead to inadequate estimations. To cope with time-variant cases, some adjustment mechanisms must be introduced in the previous basic equations. Several solutions to this problem have been proposed (Åström and Wittenmark, 1995;Ljung, 1987). When the real parameters have abrupt changes, the covariance matrix P must be periodically reset to αI, with α being a large number. If the real parameters are slowly time-variant, an RLS with exponential forgetting can be used. In the RLS with exponential forgetting algorithm, the least-squares criterion is given by, where 0 < λ < 1, designates the forgetting factor. This factor implies a time varying weighting of the data.
This method has the main disadvantage that when the input is not persistent, and as the old data is discarded in the estimation procedure, the matrix P increases exponentially with rate λ, which in this case was chosen as 0.9. This phenomenon is called estimator windup.
To solve the problem of estimator windup, an estimator that forgets the information only in the directions in which new information is gathered was implemented. For the estimator described by Eqs (9.1)-(9.8), the estimations converge avoiding large changes in the parameters. The parameters that minimise the cost function in Eq. (8) are obtained recursively from Parkum (1992): This method can be combined with other methods and techniques aiming to achieve a robust estimation algorithm. A simple way to improve the robustness of this method is to substitute for the prediction error used in the estimator law, ε(k + 1) = T i (k + 1) − ϕ T (k + 1)θ(k), with a function, f(ε(k)), in order to decrease the effects of infrequent large noise signals, such as in the case of intermittent sensor failures, over the computed parameters. In this way, the parameter adaptation law can be replaced by (Åström and Wittenmark, 1995), where f(ε(k + 1)) is a linear function for small prediction errors but increases slowly for large ones, in which d is a positive constant.

Model predictive control
Model Predictive Control (Clarke et al., 1987) comprises a collection of control methods having in common that the controller is based on the future predictions of the system behaviour using a mathematical model of the plant. There are several predictive control algorithms based on process models. These algorithms differ from each other only in the system or disturbances model's structures and in the objective function to be minimised (Camacho and Bordons, 1994), and have been applied to a wide set of control engineering  applications including greenhouse environmental control (Pinon et al., 2000;Kyriannakis et al., 2002;El Ghoumari et al., 2002).
The performance of MPC depends largely on the accuracy of the process model. This performance increases as the process-model mismatch decreases. The estimated model must be as simple as possible and capable of describing the system dynamics in a way to predict, with some precision, future outputs. So, a large part of the design effort is related to system modelling and identification. MPC involves the computation of a sequence of future control values for which it is expected that the system output tracks a given input reference. The methodology underlying these types of controllers is characterized by the strategy illustrated in Fig. 4 (Camacho and Bordons, 1994).
Future outputs for a prediction horizon L are computed for each sample k using the process model. The predicted outputŷ(k + j|k) for j = 1, . . ., L is computed using past inputs and outputs as well as the computed future values of the control signal. The collection of future control signals is computed by optimising a predefined criterion in order to maintain the process output as close as possible to the reference w(k). This criterion normally takes the form of a quadratic function of the error between the predicted output and the set point. In most cases, the control effort is included in the objective function in order to avoid abrupt changes in the control action.
Future control actions are computed optimising a specified cost function governed by the following expressions: where ε(k + j|k) is the prediction error between the future trajectory and expected output, u(k + j − 1) represents the control effort, λ 1 and λ 2 the weighing factors for each component, h max and h min represent the maximum and minimum prediction horizon and h c is the control horizon. Constants h max and h min represent the instant limits in which it is desirable that the output follows the reference.
The reference trajectory w(k + j) is sometimes different from the real reference (Clarke et al., 1987). Normally, a soft approximation from the actual value of the output towards the known reference is considered. This approach avoids abrupt changes in the control action by means of less aggressive responses. The shaped reference w(k + j) is often approximated by using a first-order lag model as described by Eqs. (14) and (15).
with α ∈ [0, 1], j = 1, 2, . . ., and r(k) denoting the real reference signal. Kennedy and Eberhart (1995) proposed the Particle Swarm Optimisation (PSO) algorithm, conceptually based on the social behaviour of groups of organisms, such as herds, schools and flocks. As an evolutionary technique, the PSO is a population-based algorithm, formed by a set of particles, representing potential solutions for a given problem. Each particle moves through a n-dimensional search space, with an associated position vector

The particle swarm optimisation algorithm
The original PSO model integrates two types of knowledge acquisition by a particle: through its own experience and from social sharing from other population members. The former was termed cognition-only model and the latter was termed social-only model (Kennedy, 1997). The behaviour of each particle is based on these two types of knowledge and their current position regarding the search. The particle behaviour is modelled by using the following equations: in which d represents the dimension index, 1 ≤ d ≤ n, p id (t) represents the best previous position of particle i in the current iteration t, p gd (t) represents the global best position in the current iteration for a pre-defined neighbourhood type, and the terms co, so, v denote cognition-only, social-only and velocity. Parameter ϕ 1 is known as the cognitive constant and ϕ 2 is known as the social constant, which represent uniformly distributed random numbers generated in a pre-defined interval. An additional parameter was incorporated into Eq. (18) by Shi and Eberhart (1999), resulting in the relation (20): in which ω represents the inertia weight. The value given to the inertia weight will affect the type of search in the following way: a large ω will direct the PSO for a global search while a small ω will direct the PSO for a local search. This parameter can vary linearly from a larger value to a smaller value in order to make the search global in the early run and local in the end of the run. Constants ϕ 1 , ϕ 2 and ω can be interpreted as the confidence that each particle has in the current position, its own experience and its neighbours experience, respectively. The neighbourhood can be of different size and topology. Each particle can take into account either: (i) The social information from a list of particles pre-defined in the beginning of the simulated evolution. The list can incorporate all the population individuals, with an individual being able to use the best solution found by every other member. This full-connected social network structure was termed Star. In other list definitions, an individual uses only k adjacent neighbours organised in a Circle or Wheel topology (Kennedy et al., 2001). (ii) The physical information which considers distance between neighbouring individuals, evaluated using some metric definition.
The velocity is limited by a maximum, V max , meaning the maximum jump that each particle can make in one iteration. The selected value for V max should not be too high to avoid oscillations, or too low to avoid search traps. The inertia weight and maximum velocity parameters selection in the PSO algorithm was studied and reported by Shi and Eberhart (1998). Each particle position should also be located within its dynamic range [X min , X max ].
The PSO algorithm has been applied successfully to solve several optimisation problems (Kennedy et al., 2001), including control tuning and optimisation (De Moura Oliveira et al., 2002;Coelho et al., 2002a,b).

Methodology of the optimisation application
The problem addressed in this report is to control the air temperature within a greenhouse using a MPC strategy. The quadratic programming (QP) problem underlying this type of controller was solved iteratively using the PSO algorithm and the results were compared with the ones obtained by using genetic algorithms (GA) and sequential quadratic programming (SQP) (Biggs, 1975) methods.
Besides the modelling technique described in Section 2.2, the MPC requires models to predict the outside climate during the prediction horizon. Several approaches were successfully used with this propose, such as the lazy men's prediction (Van Straten, 1996). Here, auto-regressive models, described by Eq. (21), are applied to predict the outside air temperature and solar radiation (Coelho et al., 2002a,b), in which A TS is a fourth-order polynomial in q −1 , y TS the time series to be modelled and ξ is the model error. To use evolutionary algorithms as design tools within the predictive control framework, it is necessary to modify them accordingly. The prediction steps are represented by population members that correspond to genes and space coordinates in the GA and PSO algorithms, respectively. Thus, control actions u to be applied to the system in a specified future time are encoded into corresponding data structures that form the population. In each generation/epoch the best two solutions found are shifted one position toward the present instant and introduced in the population of the next generation. The size of the population must be related to the size of the search space, ensuring a sufficient number of points for the evolutionary algorithm prospect. In the present case, a population of size n = 30 was selected in order to reduce the computational time. The best population/swarm solution was used as a stop criterion; i.e., the search algorithm stops if its convergence rate does not change in 30 generations/epochs.
The quadratic programming problem with linear restrictions was solved using the PSO, GA and SQP for indoor greenhouse temperature control using the described MPC strategy. Tuning parameters for both algorithms (Tables 1 and 2) are typical values used in these types of optimisation algorithms. Further study involving different parameters is underway, which will be reported in future publications, in which a comparison will be made with some results published by (Krink et al., 2001;Ursem et al., 2002;Ursem, 2003). For control purposes, the objective function to be minimized is described by: The cost function factors, λ 1 and λ 2 , define the set-point accuracy and control effort relative weights, respectively, and are set to λ 1 = 0.6 and λ 2 = 1. The performance of each optimisation algorithm was analysed in three different aspects: set-point accuracy (23), energy consumption ( (24) and (25)), for ventilation and heating signals, respectively, and the normalized computing time (26) spent by each algorithm in relation to the maximum C T = elapsed time GA,PSO,SQP max(elapsed time GA,PSO,SQP ) (26)

Results and discussion
Simulations using two data sets, for the 5th of November of 2000 and the time period from the 1st of December to the 12th December of 2000, are presented to illustrate the performances of the PSO, GA and SQP algorithms regarding set-point accuracy, energy consumption of the ventilation and heating systems, actuator effort (standard deviations of the heating and ventilation control signals, σ H and σ V ) and computation time.  Tables 3 and 4 show the results obtained with PSO, GA and SQP algorithms for the criteria defined by equations (23)-(26) using the data sets mentioned. Tables 3 and 4 simulation results, indicate that the PSO algorithm was able to reduce the set-point tracking error in approximately 30% (27%) and 81% (77%) relative to the set-point accuracies achieved with the genetic and SQP algorithms, respectively. Simultaneously, the PSO algorithm was able to decrease the variations in the heating and ventilation control signals, reducing the effort over the actuators. Regarding the energy consumption computed with the three methods, it can be concluded from Tables 3 and 4 that they are very similar for the heating and ventilation Fig. 6. Ventilation and heating control signals computed with the PSO algorithm. These signals are scaled between 0 and 1, corresponding to actuator power inputs range, from 0 to 100% of their maximum power.
actions. The major drawback of the PSO algorithm is due to the large computation run time, about 2.8 times greater than that required by the SQP method.
Due to the difficulty in visualizing the simulation curves for the largest data set, Figs. 5 and 6 show the simulated set-point tracking and actuating signals computed with the PSO algorithm only for the 5th of November, 2000. Fig. 5 indicates a close match between the system response (T i ) and set-point (w) due to the time scale used, as can been seen from the error plot clearly indicating error amplitudes comprised between ±0.2 • C. The heating and ventilation control signals are shown in Fig. 6, scaled between 0 and 1, corresponding to actuator power inputs ranging from 0 to 100% of their maximum power. It can be seen from the air temperature and the controller output curves (Figs. 5 and 6), that the PSO multivariable predictive controller has a good performance in respect to set-point accuracy, energy saving and reduction of control effort.

Conclusions and further work
The particle swarm optimisation algorithm was proposed as an alternative method to the design of a greenhouse air temperature model predictive controller subject to constraints. The controller outputs are computed to optimise the greenhouse environment future behaviour, regarding air temperature set-point tracking and control effort reduction. The prediction horizon considered was of 1 h with a 1-min sampling period. From the simulation results, we can conclude that the PSO algorithm, using a small population size, was able to significantly reduce the set-point tracking error using a lower control effort. As an evolutionary inspired optimisation algorithm, the PSO requires a larger computation run time than a classical optimisation technique, which still restricts its use in real-time applications. However, future research involving the use of micro-sized populations and parallel computational means will be directed to real-time greenhouse environmental control. Multiple tests were also conducted using data sets for other periods during the year, and the results are in agreement with those presented in this report. Moreover, the proposed techniques will be implemented and tested not only for greenhouse temperature control, but also for air humidity and CO 2 concentration control.