Branch and Bound algorithms in greenhouse climate control
Abstract
The horticultural sector has become an increasingly important sector of food production, for which greenhouse climate control plays a vital role in improving its sustainability. One of the methods to control the greenhouse climate is Model Predictive Control, which can be optimized through a branch and bound algorithm. The application of the algorithm in literature is examined and analyzed through small examples, and later extended to greenhouse climate simulation. A comparison is made of various alternative objective functions available in literature. Subsequently, a modidified version of the B&B algorithm is presented, which reduces the number of node evaluations required for optimization. Finally, three alternative algorithms are developed and compared to consider the optimization problem from a discrete to a continuous control space.
Full text
Branch and Bound algorithms in greenhouse climate control Author: Marleen Hermelink MSc Research report Universidad de M´alaga in co-operation with Wageningen University Computer Architecture and Operations Research and Logistics Grant TIN2015-66680-c2-2-R from the Spanish state in part financed by the European Regional Development Fund (ERDF) Supervisors: Dr. E.M.T. Hendrix and Dr.R. Haijema June 28, 2016
Branch and Bound Algorithms in Greenhouse Climate Control Marleen Hermelink Student ID: 92092329050 MSc program: Master in Biosystems Engineering Specialisation: Operations Research and Logistics Supervisor: Eligius Hendrix June 28 , 2016 Abstract The horticultural sector has become an increasingly important sector of food production, for which greenhouse climate control plays a vital role in improving its sustainability. One of the methods to control the greenhouse climate is Model Predictive Control, which can be optimized through a branch and bound algorithm. The application of the algorithm in literature is examined and analyzed through small examples, and later extended to greenhouse climate simulation. A comparison is made of various alternative objective functions available in literature. Subsequently, a modied version of the B&B algorithm is presented, which reduces the number of node evaluations required for optimization. Finally, three alternative algorithms are developed and compared to consider the optimization problem from a discrete to a continuous control space. Keywords: Branch and Bound, Model Predictive Control, Optimization, greenhouse climate, control action, control sequence, state variable. 1
Contents 1 Introduction 4 1.1 ResearchQuestions........................ 5 2 Literature Background 6 2.1 Climate in the Greenhouse . . . . . . . . . . . . . . . . . . . . 6 2.2 Common Control Systems . . . . . . . . . . . . . . . . . . . . 7 3 Model Predictive Control (MPC) 9 3.1 GeneralOverview......................... 9 3.2 Mathematical Notation . . . . . . . . . . . . . . . . . . . . . 9 3.3 General Problem Statement . . . . . . . . . . . . . . . . . . . 12 3.4 ExampleProblem......................... 12 3.5 SolvingMPC ........................... 14 4 Solving MPC by B&B 16 4.1 The Branch and Bound Algorithm . . . . . . . . . . . . . . . 16 4.2 General Problem Statement . . . . . . . . . . . . . . . . . . . 16 4.3 Solving by Branch and Bound . . . . . . . . . . . . . . . . . . 20 4.3.1 Dening the Search Tree . . . . . . . . . . . . . . . . . 20 4.3.2 Bounding the Search . . . . . . . . . . . . . . . . . . . 22 4.4 ExampleProblem......................... 25 5 Extension of the B&B Implementation 30 5.1 ClimateModel .......................... 30 5.2 Discretization of the Control Space . . . . . . . . . . . . . . . 32 5.3 Objective Function . . . . . . . . . . . . . . . . . . . . . . . . 34 5.3.1 Objective Function A: Sum of Squares . . . . . . . . . 34 5.3.2 Objective Function B: Sum of Costs . . . . . . . . . . 35 5.3.3 Objective Function C: Water and Energy Reduction . 37 6 Implementation of the Extensions 40 6.1 General Problem Statement . . . . . . . . . . . . . . . . . . . 40 6.2 ExampleProblem......................... 40 6.3 Results............................... 42 6.4 Design of Experiments . . . . . . . . . . . . . . . . . . . . . . 44 7 Modication for Node Reduction 47 7.1 General Approach . . . . . . . . . . . . . . . . . . . . . . . . 47 7.2 Results............................... 49 2
8 From Discrete to Continuous Control Space 52 8.1 Multi-modality of the Objective Functions on the Control Space 52 8.2 Computing the Lower Bound . . . . . . . . . . . . . . . . . . 53 8.3 Method 1: Multi-Stage Bisection . . . . . . . . . . . . . . . . 55 8.4 Method 2: Multi-stage Multisection . . . . . . . . . . . . . . 59 8.5 Method 3: Single-stage Bisection . . . . . . . . . . . . . . . . 62 8.6 Comparison of Methods . . . . . . . . . . . . . . . . . . . . . 64 9 Discussion and Conclusion 67 9.1 Discussion............................. 67 9.2 Conclusion ............................ 67 9.3 FurtherResearch......................... 69 A Appendix: Detailed Climatic Model 74 A.1 ExtendedNotation........................ 74 A.2 Complementary Climatic Model Equations . . . . . . . . . . 78 3
1 Introduction In the past two decades, horticultural production has undergone a technological revolution. Producers around the world have become increasingly competitive with help of the new available technologies, and so also have producers in the Mediterranean area [3]. Where a few years ago 100 tons of tomato per hectare was considered impressive, now a harvest of 300 tons per hectare is quite standard [3]. This exorbitant productivity has been achieved mainly through the implementation of greenhouses. These are considered ideal for growing crops, as they provide an enclosed environment which allows for controlled climate and fertigation [28]. Due to this highly productive nature, the greenhouse industry in the Mediterranean has become increasingly important to ensure global food security. Spain, for example, is often referred to as La Huerta de Europa (the kitchen garden of Europe), producing approximately 50% of all lettuce, 30% of all tomatoes, and 18% of all vegetables consumed in the EU [34]. These gures show the scale and importance of the greenhouse industry, hence clarifying the extensive research which is currently performed on it. Ongoing research on greenhouse production is very broad, but a large fraction of the issues being investigated are related to increasing the sustainability of the production system. Matters such as improved water use eciency, reduced risks of pests and diseases and reduction of chemical residues on the crop and in the soil are all becoming increasingly important issues which the industry must address. Simultaneously, they must keep costs low and product quality high [25]. Therefore, one of the main issues requiring research at the moment is energy consumption [3]. Despite the mild climate, greenhouses in the Mediterranean still require some degree of heating in the winter and cooling in the summer [9]. A reduction of energy demand for these seasons could lower the sector's production costs while simultaneously improving its environmental performance. Furthermore, due to the increasing concern on global warming, it is expected that regulations concerning CO 2 emission will become more stringent [3], so that a decrease in energy consumption will likely be necessary for the greenhouse sector. One of the main methods which has been implemented to lower energy consumption is climate control [28]. This consists of regulating the greenhouse climate so as to avoid extreme conditions which can damage the crop and to achieve suitable temperature integrals to speed up crop development and improve its quality [28]. Accurate and optimal climate control contributes to avoiding unnecessary heating and cooling, hence reducing energy consumption of the greenhouse. However, achieving this is not a simple task. Extensive research has been performed already on the control and optimization of greenhouse climate, yet many of the methods proposed thus far have some computational or practical drawbacks [7]. Therefore, it is interesting to further investigate an alternative method which has emerged during the 4
past few years, which makes use of the Branch and Bound algorithm to optimize climate regulation. This research focuses on the possibilities of this new approach, and aims to evaluate its potential as a climate control mechanism. In order to do this, the research questions in the following section are investigated. 1.1 Research Questions Main Research Question: To what extent can Branch and Bound algorithms be applied to optimize climate control in greenhouses? Specic Research Questions: 1. How has climate control in Mediterranean greenhouses been achieved thus far? 2. What kind of investigation has already been done on the topic of branch and bound algorithms for climate optimization? 3. How do the Branch and Bound algorithms for greenhouse climate optimization work? 4. What are the advantages and disadvantages of this optimization method? 5. How can the current algorithms be improved? 6. Can the algorithm be applied on a real data set from a greenhouse? 5
2 Literature Background 2.1 Climate in the Greenhouse In order to delve into the application of the Branch and Bound algorithm for climate regulation, rst some understanding of the greenhouse climate is required. The climatic conditions in the greenhouse are of great importance to the crop production. The climate aects not only the yield of the crop, but also the quality of the products [28]. To maximize the economic benet of the horticultural farm, a balance must be found between improving production and the costs of obtaining the right climatic conditions [28]. Most of the crops grown in greenhouses are adapted to temperatures between 17-27 C with a lower and upper limit of 10 C and 35 C [3]. They require a humidity within a range of 60-80% [6]. Temperatures outside this range lead to suboptimal crop production and even to permanent crop damage [1]. Too high humidity levels can lead to the development of fungi on the crop, while a humidity that is too low can cause water stress [3] [28], both of which lead to a decrease in production. Inside the greenhouse, the climatic conditions which can be controlled are temperature, humidity, Photosynthetically Active Radiation (PAR) and CO 2 concentration [28]. Temperature is the condition that inuences crop growth most directly, and is thus traditionally the main focus of climate control inside the greenhouse [3][28]. On the other hand, humidity has an indirect eect on crop growth through its inuence on crop transpiration, and should thus also be taken into consideration. However, humidity and temperature are highly inversely correlated [28]. To address this, the general solution is to keep temperature as the main control variable, but to adjust the desired temperature depending on the relative humidity [28]. The climate inside the greenhouse is aected by the outside climatic conditions (such as air temperature, wind speed, humidity, etc.) as well as by dynamic processes inside [25]. These factors are considered disturbances, as they can cause the temperature inside the greenhouse to deviate from the desired set point [28]. A complete overview of the main disturbances is given in Figure 1. To correct the temperature for these disturbances, the greenhouse has some main actuators: heating, ventilation, shading screens and fog system [28],[5]. Heating can be applied through hot air, used to avoid sudden temperature drops, or through hot water, which is applied in a more permanent way [9]. Ventilation can be performed passively, by opening windows, or forced, which is done by bringing in or extracting air with fans [9]. Screens are placed above the crop and can be folded or unfolded to create or remove shade and the fog system can be turned on to increase the humidity and reduce the temperature. Control engineering allows the automatic control of these actuators to keep track of the reference values despite the disturbances acting on the system [27]. 6
Figure 1: Inputs and outputs of climatic conditions in the greenhouse. Taken from [28]. 2.2 Common Control Systems The main goal of climate control is to keep the temperature inside the greenhouse at the desired set-point. To do this, climate control systems must determine what the best control action is at a certain moment in time. A control action can be, for example, the opening of a vent, or turning on the heating. Due to the complexity of the climate, the changes caused by actuators cannot be fully described by linear models [28]. Despite its complex nature, a number of control systems have been developed which enable climate regulation. In order to understand how climate control is generally achieved, a brief overview is given of the most commonly used systems: Proportional Integral Derivative (PID) controllers: These have been widely used for greenhouse climate control [28]. This technique is based on simplied transfer function models obtained from reaction curve tests. Three parameters must be determined beforehand (proportional gain, integral time, and derivative time) to tune the controller. In general, the goal of the controller is to quickly correct temperature deviations from the set-point. However, due to the dynamic nature of the greenhouse climate, it is dicult the obtain a good performance using xed parameter values. Therefore, this technique is often combined with other control schemes by actively adjusting its parameter values or combining it with a feed-forward loop to account for measured disturbances before they alter the climate [28]. Gain Scheduling controllers: This type of controller actively changes the control parameters based on a table or function calculated previously. 7
A drawback is that the construction of the table or function relating the parameters with the measured variables requires extensive simulation. Furthermore, there are little results on robustness, performance or stability of the system [28]. Feed Forward controllers: This type of controller is based on physical laws and measured data. The system measures the disturbances to the greenhouse climate and tries to compensate their eect before they have caused a deviation from the set-point. However, this system requires a mathematical model of the process, which is dicult to build accurately, and is thus often combined with a feedback system [28]. Besides the previously described methods, there are many other types of more complicated control systems. Some of these are model predictive control, robust control, adaptive control, optimal control and many more, as well as combinations of the dierent systems [28]. However, the focus of this report is on model predictive control, as this type of control can make use of Branch and Bound algorithms. For a detailed explanation of the other types of control, the book [28] can be consulted. 8
unsuitable for systems with short sampling times, as is the case in a greenhouse. Additionally, the convergence can lead to local minima, causing poor performance of the MPC algorithm [20]. These issues have led to the interest in using alternative methods of optimization for non-convex optimization problems [20]. These can be employed when the control space is discretized, so the problem is transformed into a discrete optimization problem. In that case, techniques such as dynamic programming [2], genetic algorithms [5, 24], and branch and bound (B&B) [4, 7, 15, 20, 30, 33] can be used. When the B&B method is applied, the discretized control space is structurally searched through a tree structure and bounds are applied to restrict the branching and avoid an enumerative search [7, 20]. This method has proven to give better results than iterative optimization techniques [33]. This is due to some of the algorithm's intrinsic properties, which lend it advantages over other techniques when applied to MPC. Previous studies point out the following advantages: The global optimum is always found, hence guaranteeing the optimality of the controller (within the discrete decision space)[33]. The method deals with constraints implicitly, and is not negatively affected by them. Constraints may even improve the bounding eciency by eliminating more branches [7]. The algorithm does not require an initial guess of the optimal solution. Hence, as opposed to iterative optimization, its performance cannot be negatively aected by poor initialization [33]. In order to test these properties and to investigate the dynamics and possible improvements, the implementation of B&B to MPC will be explored further in the next chapters. The purpose of this report is to examine the possibilities of a solution method to one time interval, regardless of where that interval is in the total time line. Therefore, from now on it is assumed that the algorithm is at a certain time interval k , but the value of k will not be specied anymore, as it is not relevant to the optimization process of a sequence V . 15
4 Solving MPC by B&B In this chapter, rst the general concept of the Branch and Bound algorithm is described in Section 4.1. Subsequently, the optimization problem is specied in more detail in Section 4.2, and Section 4.3 is an analysis of how the general concept of B&B can be applied to solve the optimization problem. Finally, the application of the algorithm is illustrated with a small example problem in the last subsection. 4.1 The Branch and Bound Algorithm The branch and bound method is characterized by the application of four rules: Branching, Bounding, Elimination and Selection [22]. In doing this, the method avoids visiting branches which are known not to be optimal [12]. The algorithm for a minimization problem starts with a set C 1 containing all feasible solutions. The Branching rule consists of splitting the set C 1 into multiple subsets. For every subset j a lower bound J L j of the minimum objective function value is determined. The lower bound is the best possible objective function value that this subset could potentially oer [12]. At every stage, there exists also a global upper bound J U of the minimum objective function value over the total feasible set. The upper bound is dened by the objective value of the best (lowest) feasible solution found thus far [12]. The Bounding and Elimination rules consist of discarding all subsets for which J L j > J U , as they can never contain a better solution than the one already found for J U . In the case that J L j < J U , then for some type of problems a feasible solution x m in the subset j can be evaluated to determine J j . If J j < J U , then J j will replace J U as the upper bound during further branching [12]. Finally, the Selection rule is to choose one of the subsets to further branch into. This inuences the performance of the algorithm, as it denes how the decision tree will be searched [12]. For example, choosing the subset with the lowest lower bound for further branching yields a depth rst search, while choosing the largest subset yields a breadth rst search. The main target is to quickly nd a sharp upper bound J U , so that large parts of the tree can be pruned [12]. 4.2 General Problem Statement Although the application of the B&B algorithm to MPC does follow the same rules of Branching, Bounding, Elimination and Selection normally used in B&B optimization, it's search structure diers from the traditional algorithm described in the previous section. In order to solve control sequence V with branch and bound, it is necessary to dene the objective function J of the control problem cumulatively over the prediction horizon [5, 7, 20]. Each control action v i of a control sequence V contributes a certain transi16
tion cost c to its total objective cost J ( V ): J ( V ) = PH X i =1 c ( v i ; s i ) with V = ( v 1 ;:::; v PH ) (9) The dependence of the objective value on V is henceforth considered implicit in the simplied notation J = J ( V ) Up to the i th control action of V , the partial cost thus far can be calculated by adding the transition costs up to i : J i = i X ` =1 c ( v ` ; s ` ) 8 i = 1 ;:::;PH (10) Therefore, the cost J PH up to control action v PH is equal to the total costs J of the whole control sequence: J PH = PH X ` =1 c ( v ` ; s ` ) J PH = J (11) The transition cost c ( v ; s ) is a function dependent on the control action v and on the state s at that interval. It computes the cost of implementing a certain control action v , as the actuator settings. The state vector s i contains the state of the system at the beginning of interval i . It contains NS elements, consisting of two dierent types of state variables. The rst elements are climatic state variables, which are states of the climate that must be controlled or followed. Examples are temperature, humidity, ground temperature, or the PAR radiation in the greenhouse. The number of climatic state variables is denoted NSC . The last elements of s i are actuator state variables, as they contain the control action v i 1 that was implemented in the previous interval. This works as the memory of the algorithm, in order to be able to recall the previous control action. Since v i 1 contains NA elements, the last NA elements of s i together comprise v i 1 . The total number of elements NS in the state vector s i is thus NS = NSC + NA The state vector s i is hence dened as: 17
s i = 0 B B B B B B B B @ s 1 . . . s ( NS NA ) s ( NS NA +1) . . . s NS 1 C C C C C C C C A = 0 B B B B B B B B @ s 1 . . . s ( NSC ) s ( NSC +1) . . . s NS 1 C C C C C C C C A 9 > = > ; Climatic state variables 9 > = > ; v i 1 (12) 8 i = 1 ;:::;PH The transition cost function c ( v ; s ) is a topic on which the available literature diers, so this will be expanded on in depth in Section 5.3. For now, a conventional MPC objective function will be used, which uses sumquadratic functions to minimize the overshoot (error) and the control eort [32]. The error e must be minimized to maintain the necessary conditions to insure the quality of the products [28]. On the other hand, the control eort, denoted v , is minimized, because every change in the control action requires energy [5]. The transition cost is calculated as follows: c ( v ; s ) = jj e jj 2 1 + jj v jj 2 2 (13) The error e is dened as the dierence between the setpoint SP and the predicted state s i +1 = T ( s;v ) at the beginning of the next interval [20], where T is the transition function of the state: e = T ( s;v ) SP (14) The setpoint vector SP contains the desired values of the climatic state variables of the greenhouse. Its last elements have a value of zero because SP must have the same dimensions as the state vector s to compute Eq.(14). Since s ( NSC +1) to s NS are not related to the climatic state of the greenhouse, they do not have a setpoint, hence leading to zero elements in SP : SP = 0 B B B B B B B B @ SP 1 . . . SP ( NSC ) SP ( NSC +1) . . . SP NS 1 C C C C C C C C A 9 > = > ; Climatic state setpoints 9 > = > ; = 0 (15) The control eort v at interval i is the change in control action from the previous interval i 1 to i , with i = 1 ;:::;PH . Recall from Eq.(12) 18
that the control action at i 1 are stored as the last elements s ( NSC +1) to s NS in the state vector s . This functions as the "memory" of the algorithm to recall the previous control action v i 1 . Therefore, the control eort v is dened as: v i = v i v i 1 8 i = 1 ;:::;PH (16) v = v ( s ( NSC +1) ;:::;s NS ) T (17) Finally, 1 and 2 are matrices which simultaneously normalize and weigh the error and control eort [20]. For the normalization, the matrices must map the terms linearly in the range [0,1] by dividing each element by its user-estimated maximum. The weights W 1 and W 2 of the terms should be a value between 0 and 1, and together should sum to 1 [29]. The weights are used-dened, and are meant to reect the priorities or goals of the user. The matrices look as follows: 1 = W 1 0 B B B B B B B B @ e 1 1 0 0 0 . . ..... . . 0 . . . 0 e NSC 1 0 0 0 0 0 0 0 . . .. . .. . .. . ..... . . 0 0 0 0 1 C C C C C C C C A (18) 2 = W 2 0 B B B B B @ v 1 1 0 0 0 v 2 1 0 . . .. . ..... . . 0 0 v NA 1 1 C C C C C A (19) where e 1 ;:::;e NSC are the maximum errors of climatic state variables, and v a the maximum change of control variable v a , with a = 1 ;:::;NA . Note that the last rows of 1 contain only zeros because the last elements of s comprise the control action at the previous interval i . Since these are climatic state variables, they are not relevant to the deviation term e . The general optimization problem statement for a control sequence V can thus be formulated as follows: 19
min V J ( V ) = PH X i =1 c ( v i ; s i ) with V = ( v 1 ;:::; v PH ) v = ( v 1 ;:::;v NA ) s = ( s 1 ;:::;s NS ) c ( v ; s ) = jj T ( v ; s ) SP jj 2 1 + jj v s ( NSC +1 ;:::;NS ) jj 2 2 s.t. s i +1 = T ( v i ; s i ) 8 i = 1 ;:::;PH v i = v CH 8 i = CH + 1 ;:::;PH v ia 2 a 8 i = 1 ;:::;PH 8 a = 1 ;:::;NA 4.3 Solving by Branch and Bound The optimization problem outlined in the previous section can be solved using a branch and bound search method. The search method is analyzed in the following subsections. 4.3.1 Dening the Search Tree To implement the B&B to solve the optimization problem, rst the search tree must be dened. Each time interval i , with i = 1 ;:::;PH , represents a level in the search tree ( i = 0 at the initial node) [20]. At each level i within the control horizon CH , i = 1 ;:::;CH , the algorithm must decide the control action v i . The system has B alternatives for v i . Therefore, each node will have B branches [20]. The j th branch, with j = 1 ;:::;B , is represented by ! j . For example, consider a greenhouse with a heater as single actuator v 1 , which can be turned o or on, so v 1 2 f 0 ; 1 g . Since the problem has only two control alternatives, at each level i of the tree each node will branch in two. So v i can be either 0 or 1. Therefore, the algorithm can choose between ! 1 and ! 2 for v i , where control alternative ! 1 = (0) is to turn the heater o and control alternative ! 2 = (1) is to turn the heater on [20]. The transition cost c of implementing branch ! j as control action v i is c ( ! j ; s i ), calculated by Eq.(13). Although for this simple problem this notation might seem overly complicated or redundant, more complex problems later on will justify its necessity. Figure 4 gives a graphical representation of a search tree with two control alternatives. Figure 4 also shows that no branching takes place beyond the control horizon ( i > CH ). Control action v CH is applied successively until PH . In other words, when level CH is reached, the control action at the last interval within the control horizon is applied [20]. So far the tree has been illustrated in the case of only one actuator v 1 . However, in reality a greenhouse has more than one actuator. In that case, 20
Figure 4: Example of a branch and bound tree applied to a MPC optimization problem with one actuator v 1 2 f 0 ; 1 g , prediction horizon PH = 5, and control horizon CH = 3. the tree will look largely the same as for the case of one actuator, except that the number of control alternatives B will be much larger. Recall from Eq.(4) that the control action v i contains NA control variables, one for each actuator. Without loss of generality, it is assumed that each actuator v a with a = 1 ;:::;NA , has M discrete alternative control actions. Therefore, a discrete control action b of actuator v a is represented by ! ab [20]. The set of all possible discrete control actions for the actuator v a is represented by a [20]: a = f ! ab j b = 1 ;:::;M g (20) The discrete set containing all the possible control actions is hence given by [20]: = 1 2 NA (21) The number of the total possible discrete control alternatives B is given by [20]: B = M NA (22) Each vector of the matrix can be represented by ! j , with j = 1 ;:::;B [20]: = f ! 1 ;! 2 ; ;! B g (23) Each vector ! j is therefore a control alternative of the actuators, and corresponds to the j th branch at each tree node. Hence, at each time interval, 21
B control alternatives can be implemented for v i , resulting in a maximum of B branches per node [20]. The general problem is depicted in Figure 5. Figure 5: General case of a branch and bound tree applied to MPC, with a prediction horizon PH , control horizon CH , and B control alternatives for control actions v i with i = 1 ;:::;PH . Taken and modied from [20]. As an example, consider a greenhouse with actuators v 1 and v 2 , both of which have M = 3 control alternatives. Assume that v 1 2 f 0 ; 50 ; 100 g and v 2 2 f 0 ; 30 ; 70 g . The maximum number of branches will be B = M M = 3 3 = 9. The problem will thus have B = 9 branches or control alternatives per node. The matrix of control alternatives is given as an illustrative example in Table 2. Table 2: Control alternatives ! 1 ;:::;! 9 for a problem with actuators v 1 2 f 0 ; 50 ; 100 g and v 2 2 f 0 ; 30 ; 70 g . ! 1 ! 2 ! 3 ! 9 v 1 0 0 50 100 v 2 0 30 70 70 4.3.2 Bounding the Search The characteristic of the branch and bound algorithm that distinguishes it from an enumerative search is the Bounding rule. If all the possible branches of the tree were searched, it would result in B CH dierent control sequences, which even for a small number of actuators and discretizations 22
can be too large [20]. Hence, the bounding is essential for the applicability of the algorithm. To avoid extensive data storage, the algorithm performs a depth rst search. Therefore, a branch ! j at level i is evaluated to determine if it can be followed to the next level i + 1 of the tree. The Bounding rule states that a particular branch is only followed if the previous cumulative cost J i 1 plus a lower bound on the cost from the level i to PH , denoted J ( i ) L , is lower than an upper bound of the total cost, denoted J U [20]. In other words, at each branch, a lower bound J ( i ) L must be computed of the estimated costs of the current and the remaining intervals, see Figure 6. The lower bound is assumed to consist of the transition cost c ( ! j ; s i ) at interval i plus the cost of the remaining intervals from i +1 to PH , see Figure 6. However, the cost of the remaining intervals is very dicult to estimate. As the lower bound may only be an under-estimate, the remaining cost is set to 0, to avoid bounding branches which could contain a better solution [7, 20]. Therefore, the bounding condition is the following [20]: J i 1 + J ( i ) L < J U J i 1 + c ( ! j ; s i )+0 < J U J i < J U (24) When the terminal level of the tree i = PH is reached, the cumulative cost J i is the terminal cost J PH . At the bottom of the tree, the path followed down the tree is a complete control sequence V . The objective value J ( V ) of the sequence is J ( V ) = J PH = J (see Eq.(11)). If this terminal cost J of the sequence is lower than the current upper bound J U , then J U is replaced by J and the control sequence V is saved as the best so far. After the whole tree has been evaluated, the control sequence stored at that moment as best so far can be concluded to be the overall optimal sequence. The initial value of J U can be set arbitrarily. In theory, it could be set to innity, as it will be replaced by a better value by the rst branch reaching the terminal level of the tree. However, to decrease the number of computations, the upper bound should always be as low as possible. Therefore, nding an initial upper bound close to the optimum can be convenient. In most of the literature, a greedy algorithm is used for this [7, 20]. The rst path followed through the search tree chooses the smallest transition cost c ( ! j ; s i ) at each level i . This strategy yields a terminal cost J PH close to the optimum and can hence reduce the necessary iterations to nd the optimal solution. The logical steps that the algorithm takes to evaluate the search tree are shown in Algorithms 1 and 2. 23
Figure 6: Graphical representation of upper bound J U , cumulative cost J i 1 and transition cost c ( ! j ; s i ) for a problem at i = 3. Algorithm 1 ControlSequence( PH;CH;B; s 1 ) Require: PH; CH; B , and s 1 Compute alternative settings := f ! 1 ;:::;! B g f V;J U g := UpperBound ( PH;CH;B; s 1 ; ) Let a := for a = 1 ;:::;CH i := 1 and J 0 := 0 while ( a 6 = ? for a = 1 ;:::;PH ) remove ! from i calculate transition cost c ( !; s i ) and new state s i +1 calculate total cost up to i : J i := J i 1 + c ( !; s i ) if ( J i < J U ) v i := ! if ( i < PH ), i := i + 1 else update upper bound J U := J i and V := ( v 1 ;:::; v PH ) i := i 1 while ( i = ? ) i := i := i 1 endwhile if ( i > CH ), a = f v CH g for a = i;:::;PH endwhile return: V , J U 24
s i = 0 B B B B B B @ s 1 s 2 s 3 s 4 s 5 s 6 1 C C C C C C A ! Air Temperature ! Air Relative Humidity ! Ground Temperature ) v i 1 (26) The climate model was taken from articles [5] and [13]. In state space form, the model is dened as a mass and energy balance [5]. The state space equations of the climate model are dened in Eq.(28) to (29). The symbolic representation of the equations was adapted to show clearly which are the dependent and independent variables. All capital letters are dependent variables, all symbols represent constant parameters, and non-capital letters are the independent variables. Since all the parameter values of the equations are xed, the model can be said to be deterministic. The actuator and state variables, as well as the disturbances which the model takes into account are given in Table 4. d s 1 d t = Q 2 ( d 3 ) Q 3 ( s 1 ;d 1 ) + Q 7 ( s 1 ;s 3 ) F 7 ( s 1 ;s 2 )( Q 4 ( s 1 ;s 2 ) + Q 6 ( v 2 ;s 1 )) Q 5 ( v 1 ;s 1 ;d 1 ) + Q 1 ( v 3 ) = ( 18 24 5 ) (27) d s 2 d t = f 3 ( s 1 ;G 2 ( v 1 ;s 1 ;s 2 ;d 2 ) + F 7 ( s 1 ;s 2 )( G 1 ( v 2 ) + G 3 ( s 1 ;s 2 )) 24 18 ) (28) d s 3 d t = Q 8 ( d 3 ) Q 7 ( s 1 ;s 3 ) Q 9 ( s 3 ) 2 4 (29) The data of the model has a sample time of 15 seconds [5]. Consequently, for the optimization problem, time was split into time intervals of 15 seconds. Within that time interval, it is assumed that the change in state is linear. Therefore, the changes in state over one interval i are: s 1 = 15d s 1 d t s 2 = 15d s 2 d t s 3 = 15d s 3 d t The transition function T of the problem is therefore dened as: 31
T ( v ; s ) = 0 B B B B B B @ s 1 + s 1 s 2 + s 2 s 3 + s 3 v 1 v 2 v 3 1 C C C C C C A (30) Table 4: Climatic state variables, disturbances and control variables (actuators) taken into account in the climate model of [5] and [13] . Climatic State Variables s 1 C Inside temperature s 2 % Inside relative humidity s 3 C Ground temperature Disturbances d 1 C Outside temperature d 2 % Outside relative humidity d 3 Wm 2 Outside Solar radiation d 4 ms 1 Outside Wind speed Control Variables v 1 % Window opening v 2 % Fog system v 3 % Heating The climate model was developed with data from a rose hydroponic crop in a plastic greenhouse with arch shaped roofs located in Moncada (Valencia, Spain) in the summer [5]. Thus, the model parameters are specic for these greenhouse conditions. A detailed account of the used parameters as well as the complementary equations can be found in Appendix A. For information on how the parameters were obtained, article [5] can be consulted. The graphs in Figure 9 illustrate the eect of each actuator on the climate state. These are meant as a visual aid to grasp the eect of the transition function. Additionally, the graphs were computed to check the quality of the climate model, as the results they show should not be counter-intuitive (for example, declining temperature with increasing heating). 5.2 Discretization of the Control Space As discussed in Section 3.5, the application of branch and bound to a control sequence optimization problem requires the discretization of the control space. Hence, an actuator v 1 such as a heater for example, which has a con32
(a) Temperature (b) Humidity Figure 9: Temperature and humidity after implementation of the range of control actions during one time interval, assuming parameters and outside weather conditions given in Table 11 in the Appendix tinuous range between 0 and 100%, could be discretized to v 1 2 f 0 ; 50 ; 100 g , implying that the heater can only be turned o, half power, or full power [5]. The drawback of the discretization is that it can cause the output to oscillate around a reference trajectory, instead of following it exactly [33]. However, literature shows that crops are more responsive to the average daily temperature than to accurate temperature evolution during the day [14, 35]. This led to the concept of Temperature Integration (TI) [17], where the greenhouse temperature is allowed to uctuate within certain boundaries, as long as the average over a certain period is maintained. Hence, the application of TI to a greenhouse allows temperature oscillations caused by the control space discretization. Recall from the previous section that the actuators that were implemented in the programmed model are the window opening v 1 , the fog system v 2 , and the heating v 2 . For the discretization, each actuator's physical limitations of saturation (maximum) and output resolution (minimum step size) were also considered. For example, a window can only be opened between 0 and 100% (saturation) and the teeth of the window rack allow a minimal movement of 5% (output resolution)[28]. Therefore, all the elements in the control space of the window opening should be within 0 and 100% and must be represented by multiples of ve. Another restriction on the discretization was the computation time, since a slight increase in the number of alternative control actions causes a drastic increase in search branches and thus also in computational time. Therefore, each actuator was restricted to a maximum of four discrete control settings. The settings were spread evenly across the total range of 0 to 100%. Finally, 33
as Section 5.1 indicates, the climate model has parameters for the summer period, so it is assumed the heater is turned o. Actuator v 3 was thus restricted to 0. The actuators of the model were discretized as follows: v 1 2 f 0 ; 33 ; 66 ; 100 g v 2 2 f 0 ; 33 ; 66 ; 100 g v 3 2 f 0 g 5.3 Objective Function In order for the optimization problem to be relevant, its objective function J should accurately reect the goals of the user. Since the literature on this topic diers, the issue was further investigated by implementing three dierent objective functions chosen from relevant articles. All three functions work according to the same cumulative principle outlined in Eq.(9), in Section 4.2, where the transition cost c at each interval i , with i = :::;PH , is added to form the total cost J of a control sequence V . Moreover, all three functions compute the transition costs c as a summation of normalized and weighed terms, each term representing an objective of the user. However, the content of the terms diers per objective function, depending on the overall approach of the function. The dierent objective functions are described in the next sub-sections. 5.3.1 Objective Function A: Sum of Squares The rst objective function is the same that was used for the short example in the previous section. It is considered the classic MPC objective function and was taken from article [20]. The objective function seeks to minimize the deviation, or error, e from the setpoint SP and the control eort v . This function was described in Section 4.2, so this chapter can be referenced for details. As a brief refresher, the equation of the transition cost is described below: c ( v ; s ) = jj e jj 2 1 + jj v jj 2 2 with e = T ( v ; s ) SP and v = v ( s 4 ;s 5 ;s 6 ) T The normalization and weight matrices 1 and 2 are dened according to the following equation, with e 1 and e 2 the maximum temperature and humidity deviation from the setpoint, and v 1 ; v 2 and v 3 the maximum control eorts of the window, fog system and heater respectively. It is 34
assumed that the minimum change of the control action is always 0. 1 = W 1 0 B B B B B @ e 1 1 0 0 0 0 e 2 1 0 0 0 0 0 0 . . .. . .. . ..... . . 0 0 0 0 1 C C C C C A 2 = W 2 0 B B @ v 1 1 0 0 0 v 2 1 0 0 0 v 3 1 1 C C A (31) Notice that the last rows of 1 contain only zero elements. The third row is empty because the third state variable s 3 is the ground temperature, which is not a state control variable. It is therefore not relevant to the deviation from the setpoint, since the user does not aim to control it. All other rows are empty because the state variable elements s 4 to s 6 are the previous control action v i 1 , which is also irrelevant to the deviation. Recall that all three actuators have the same range of 0 to 100% (see Section 5.2). Therefore, v 1 = v 2 = v 3 = 100 (32) Consequently, 2 can be simplied to: 2 = W 2 100 An important characteristic of this objective function is that its deviation term e includes both the temperature and the humidity. The main advantage of this feature is that it allows the user direct control over both the temperature and the humidity inside the greenhouse, whereas normally only the temperature is controlled. However, a disadvantage of the function is that it is not concerned with objectives such as water or energy use minimization, even though these are also important to achieve a more sustainable production. 5.3.2 Objective Function B: Sum of Costs The second objective function was taken and adapted from [7]. Similarly to function A, it strives to minimize the deviation from the setpoint SP as well as the control eort v . The function also minimizes the energy and water use. Furthermore, besides having a setpoint, the function also has a range within which the states are allowed to uctuate. Any solution going outside the range receives a penalty. The transition cost c for each level of 35
the search tree is calculated as follows: c ( v ; s ) = 1 e + Error + 2 NA X a =1 v a + ) Control Eort + 3 NA X a =1 v a w a + ) Water and Energy Use + 4 P ( v ; s ) Penalty (33) The rst term of the equation is the cost of the deviation of the states with respect to the setpoint. However, this objective function considers only the temperature as a climate state control variable. The states of the humidity and ground temperature have no setpoint, and hence do not contribute to the objective value. As a consequence, their corresponding rows in the normalization and weight matrix 1 only contain zero valued elements. The matrix is dened as: 1 = W 1 0 B B B @ e 1 1 0 0 0 0 0 . . .. . ..... . . 0 0 0 1 C C C A The second term of the objective function is the costs of the control eort, similarly to function A. However, in this case the changes in actuator settings are all added and then normalized and weighed by 2 , which is a single value: 2 = W 2 v 1 + v 2 + v 3 1 Making use of Eq.(32), 2 is simplied to: 2 = W 2 300 (34) The third term is meant to minimize water and energy consumption by attributing a cost to the use of each actuator. The importance attributed by the user to each actuator is indicated by the weights in vector w . For example, as opening a window requires much less energy and water than the heater and fog system, it is expected that the lowest weight is attributed to this actuator. All the elements in w should be in the range [0,1], and together should sum to 1. Weight vector w is dened as: w = 0 @ w 1 w 2 w 3 1 A (35) 36
where w 1 ;w 2 and w 3 are the weights of the window, fog system and heating respectively. Since all three actuators have a maximum value of 100, and the weights must sum to one, the weighted sum can never be greater than 100. Therefore, the weight and normalization value 3 is: 3 = W 3 100 (36) Finally, the last term is a penalty for the violation of the temperature and humidity ranges in which the states are allowed to uctuate (see Section 5.2). The penalty P is given by: P ( v ; s ) = max f 0 ; R T ( v ; s ) ; T ( v ; s ) R g (37) The vectors R and R are the lower and upper bound respectively of the allowed range of the temperature and the humidity. The vectors are dened as follows: R = 0 B B B B B @ T H 0 . . . 0 1 C C C C C A R = 0 B B B B B @ T H 0 . . . 0 1 C C C C C A (38) Recall that elements s 3 to s 6 of the state vector s are the ground temperature and the previous control action. As both are not relevant to the penalty term, the corresponding elements of the vectors R and R are irrelevant, and are thus set to zero. Similarly, the same rows of the matrix 4 consist only of zero-elements as well. The matrix is dened as: 4 = W 4 0 B B B B B B @ P 1 1 0 0 0 0 P 2 1 0 0 0 0 0 0 . . .. . .. . ..... . . 0 0 0 0 1 C C C C C C A (39) where P 1 and P 2 are the maximum temperature and humidity violations of the range, and are dened by the user. 5.3.3 Objective Function C: Water and Energy Reduction The third objective function's main goal is to reduce the water and energy consumption as much as possible, given that the temperature and humidity stay within a certain range. The function was taken and adapted from article [5]. The transition cost is calculated as follows: 37
c ( v ; s ) = 1 v 1 + Control Eort for Windows + 2 v 2 + Water Use (fog) + 3 v 3 + Energy Use (heating) + 4 P ( v ; s ) + Penalty + oset( v ; s ) (40) The rst term of the equation evaluates the control eort. However, only the control eort of the windows is considered, as this is the only actuator which requires a signicant energy consumption to physically change its settings [5]. The second term evaluates the water use by the fog system, and the third one considers the energy consumption by the heating system. As all three terms consider only one actuator, their weight and normalization matrices 1 , 2 , and 3 are single values: 1 = W 1 100 2 = W 2 100 3 = W 3 100 (41) The fourth term is a penalty for the deviation from the established temperature and humidity range. The penalty P is calculated using Eq.(37), dened in the previous objective function. However, objective function B allows the user to determine the normalization of the penalties, whereas [5] denes the normalization and weight with 4 : 4 = 0 B B B B B @ 15 1 0 0 0 0 20 1 0 0 0 0 0 0 . . .. . .. . ..... . . 0 0 0 0 1 C C C C C A (42) Note that a higher importance is given to the temperature (1 = 15) than to the humidity (1 = 20). Finally, the last term is an oset. This term is added to the equation to ensure that every solution going outside the allowed range has a greater objective value than those solutions inside the range. In other words, by using the oset, solutions with a penalty will always be less favorable than solutions without a penalty. The term is calculated as follows: oset( v ; s ) = ( 3 ; P ( v ; s ) 6 = 0 0 ; P ( v ; s )=0 (43) , where penatly P is given by (37). The oset has a value of 3 because the weights W 1 , W 2 and W 3 are all in the range [0,1], so the sum of the rst three terms can never be greater than 3. Therefore, the oset ensures that solutions with any violation of the ranges always have a higher value than those solutions with no violation. 38
The main advantage of this objective function is that its oset term prevents it from exploring solutions that are not promising, thus cutting o more branches of the search tree from the beginning. This should drastically reduce the computational load of the algorithm. However, as it does not have a term to follow a setpoint, it can prove to be less accurate than the other two objective functions. 39
6 Implementation of the Extensions 6.1 General Problem Statement The extensions of Section 5 were incorporated in the optimization problem modeled in Matlab. This led to an optimization problem with a more realistic state transition function, a clearly dened discrete control space, and three alternative objective functions. The general problem formulation is given below. min V J ( V ) = PH X i =1 c ( v i ; s i ) with V = ( v 1 ;:::; v PH ) v = ( v 1 ;v 2 ;v 3 ) s = ( s 1 ;:::;s 6 ) c ( v ; s ) = Eq.(13), (33) or (40) s.t. s i +1 = T ( v i ; s i ) 8 i = 1 ;:::;PH v i = v CH 8 i = CH + 1 ;:::;PH v i 1 ;v i 2 2 f 0 ; 33 ; 66 ; 100 g ; v i 3 2 f 0 g 8 i = 1 ;:::;PH 6.2 Example Problem The model built in Matlab R2015b solves the optimization problem outlined above by applying the branch and bound search structure from Section 4. The three alternative objective functions were built as separate functions, so that the model could optimize either of the three, depending on which one is called on. Since the transition function of the problem is deterministic, the search tree of a given set of weights and parameters will always result in the same solution. Indeed, because of the absence of stochasticity in the problem, it can even be said that the search tree is known beforehand. However, applying the branch and bound method avoids the need for a full enumeration of the whole tree. The model was used to solve an example problem with a prediction and control horizon of PH = 9 and CH = 5 respectively. The parameters and weights of the objective functions established for the problem are given in Tables 5, 6, 7, and 8. The parameters in Table 5 regarding the state of the climate inside the greenhouse were partially taken from literature and partially determined from personal judgement. The outside climatic conditions given in Table 6 were taken from climatic data recorded in Almera. Notice from the tables that for this particular problem, the initial temperature is outside of the permitted range, while the humidity is inside. Additionally, the outside temperature is lower than the desired setpoint, while the outside humidity is higher. 40
7 Modication for Node Reduction 7.1 General Approach The computational time required to solve the MPC optimization problem with the B&B search algorithm is caused by the large number of nodes that must be evaluated. The nodes after the control horizon are relatively simple. Those nodes have only one branch (see Figure 5), since the control action remains constant from there until the end of the prediction horizon. Therefore, this portion of the search tree could be potentially interesting for computational time reduction. This possibility was investigated creating a modied algorithm which searches the nal nodes slightly dierently. Since all the nodes of a vertical branch after the control horizon have the same control action v CH , the cost and state at each node can be simulated until the prediction horizon in one single step. Consequently, after the control horizon, the modied algorithm does not create a new node at each level i , but rather evaluates all intervals after control horizon CH in one single node. A graphical representation is given in Figure 14. The general idea behind the modication is that a reduction in number of node evaluations will lead to a reduction in computational time. To avoid confusion, the new algorithm uses a new index m for the levels of the tree instead of the index i , which has been used thus far. Up to the control horizon, both indices can be considered as alias indices, since m = i . However, in the modied algorithm, the nodes i = CH + 1 ;:::;PH are comprised into one node m = CH + 1. Therefore, m is dened as m = 1 ;:::;CH + 1. In the case that CH = PH , there are no nodes after CH , so m = 1 ;:::;CH . In general form, m can thus be dened as: m = 1 ; 2 ;:::;CH + h h = ( 0 ; if CH = PH 1 ; if CH < PH The transition cost at each node is re-dened as C . At the nodes where m 6 CH the transition cost is calculated in the same manner as has been done thus far: C = c ( ! j ; s i ), using Eq.(40). At m = CH + 1, the transition cost C is the sum of all the individual transition costs c ( v CH ; s i ) for i = CH + 1 ;:::;PH (recall that for i > CH the control action v remains constant at v i = v CH ). Therefore, the transition cost C at each level m is calculated as follows: 47
C ( !; s m ) = 8 > > < > > : c ( !; s m ) m 6 CH PH X i = CH +1 c ( v CH ; s i ) m = CH + 1 (44) with s i +1 = T ( v CH ; s i ) 8 i = CH;:::;PH The logical steps of the modied algorithm are outlined in Algorithm 3. Figure 14: Graphical representation of the modied algorithm, in which the search tree simulates the nodes after the control horizon in one single node. 48
Algorithm 3 ControlSequenceModied( PH;CH;B; s 1 )) Require: PH; CH; B , and s 1 Compute alternative settings := f ! 1 ;:::;! B g f V;J U g := UpperBound ( PH;CH;B; s 1 ; ) Let a := for a = 1 ;:::;CH m := 1 and J 0 := 0 while ( a 6 = ? for a = 1 ;:::;CH + 1) remove ! from m calculate transition cost C ( !; s m ) with Eq.(44) calculate total cost up to m : J m := J m 1 + C ( !; s m ) if ( J m < J U ) if ( m 6 CH ) v m := ! m := m + 1 else v i := ! , for i = CH + 1 ;:::;PH update upper bound J U := J m and V := ( v 1 ;:::; v PH ) m := CH while ( m = ? ) m := m := m 1 endwhile if ( m > CH ), m := f v CH g endwhile return: V , J U 7.2 Results The modied algorithm was created under the premise that a reduction in node evaluations will also lead to a reduction in computational time. To test this, both the original and the modied algorithms were used to solve the same optimization problem, and the number of node evaluations and computational time of both were compared. The control problem on which the models were tested is formulated in Section 6.2. Due to the practical issue of time availability for this report, the modied algorithm was implemented with objective function C, as this was deemed the easiest objective function in Section 6.4. To obtain a more comprehensive overview of the eect of the modication in the algorithm, the problem was also solved at dierent prediction horizon lengths. At all the tested PH values, the CH : PH ratio was kept constant at 2:3. These results are plotted in Figure 15. Additionally, the models were tested on the same problem with a constant prediction horizon of PH = 20 49
but with dierent ratios of control to prediction horizon. The results are given in Figure 16. (a) Number of node evaluations (b) Computational time Figure 15: Node evaluations and computational time of the original and the modied algorithm at increasing PH values under a constant CH : PH ratio of 2:3. Figure 15a shows that the problem evaluates a constant number of nodes at prediction horizon values up to PH = 19. This is likely caused by the fact that at these prediction horizon lengths, the algorithms can stay within the admitted temperature and humidity range without implementing any actuator change (see example in Figures 11 and 12). For a larger prediction horizon PH , the problem becomes more complex and both the number of node evaluations and the computational time increases exponentially. The graphs show that by implementing the modied model, the number of node evaluations is not reduced. Although the data did show a slight reduction, it was so small that it does not show in the graph. On the other hand, Figure 15b shows that at PH values where the number of node evaluations was high ( PH > 19), there is a substantial reduction in computational time. Figure 16 shows the result of investigating the length of the control horizon with respect to the prediction horizon. It was expected that for a longer control horizon, the complexity of the problem would be higher and thus the number of node evaluations and computational time would increase. However, both graphs show that there is no such pattern in the observed data. There is a clear peak at CH : PH = 5 : 20 = 0 : 25 and otherwise the nodes and the time stay within a relatively small range from each other. It was also unexpected that at the highest CH : PH ratios of 0.9 and 1, the computational eort was the lowest. The data shows that implementing the modied model in most cases did not lead to a great 50
(a) Number of node evaluations (b) Computational time Figure 16: Node evaluations and computational time of the original and the modied algorithms at increasing CH : PH ratio under a constant prediction horizon of PH = 20 reduction of computational eort, except at the peak around CH : PH = 0 : 2. At those points, the data shows a reasonable reduction in number of node evaluations as well as computational time. Therefore, the general overall conclusion from Figures 15 and 16 is that the adjustment to the model does lead to a computational time reduction, but only in the cases where the computational eort is high. 51
8 From Discrete to Continuous Control Space Thus far the issue of greenhouse model predictive control has been explored with a Branch and Bound algorithm applied to problems with discrete control space. However, as in reality the control space of greenhouse actuators is continuous, it is interesting to examine the possibility of applying the branch and bound algorithm to a continuous control space. This would eliminate the previously discussed need for a trade o between discretization and control performance (Section 6.4). This chapter outlines how the transition from discrete to continuous control space was approached. 8.1 Multi-modality of the Objective Functions on the Control Space In order to apply the branch and bound algorithm to a problem with continuous search space, rstly the multi-modality of each of the three objective functions of Section 5.3 must be examined. This is necessary to determine if there are multiple local minima; or in other words, if the functions are nonconvex. However, the control problem has a very high dimensionality, since the setting of each actuator at each time interval is a variable. A simple problem with, for example, a control horizon CH = 4 and three actuators, already has n = 4 3 = 12 dimensions. Since such a high-dimensional control space cannot be examined graphically, an alternative method was applied to determine the presence of multiple optima. The control space was explored by applying Multistart. This is an algorithm where local searches are performed from randomly generated starting points [12]. Each starting point will reach one of the local optimum points in the search space. The local minimum that is reached depends on the region of attraction in which the starting point is situated, and on the local optimizer that is used [12]. For all three objective functions, 100 starting points were generated and solved locally with the FMINCON function in Matlab R2015b. This was done repeatedly for increasing prediction horizons, always given a ratio of control to prediction horizon of CH : PH = 2 : 3. The resulting solutions, i.e. local minima, were rounded o to integer numbers. This implies that all control settings within a distance of 1 from each other are considered to have the same value. Although this is a relatively rough method of examination, it is sucient for the purpose of determining the presence of multiple minima. The number of local minima that resulted for each run is shown in Figure 17. The graph shows that objective function A always has only one solution, regardless of the problem's prediction horizon. The minimum point found by the local solver is thus the global optimum. Admittedly, there is a possibility that if the number of random starting points were higher, the local solver would nd more than one minimum. However, for simplication purposes 52
Figure 17: Number of local minima for objective functions A, B and C at increasing prediction horizon length using 100 random starting points and a ratio of PH : CH = 2 : 3. it will be assumed that the function has a convex shape, and thus does not require any non-convex optimization method. Function A is thus dismissed, as it is of no added value to apply the B&B algorithm to its search space. On the other hand, objective functions B and C show multiple local minima. Function B seems to be convex for problems with PH < 6, but as the dimensionality of the problem becomes higher, the number of local minima increases drastically. Conversely, function C has a large number of minima regardless of the length of the prediction horizon. With increasing prediction horizon PH , the number of minima increases slightly, yet not as steeply as for function B. The histograms in Figure 18 show that for both functions B and C, the local minima are distributed over dierent objective function values, hence proving that they are truly dierent from each other. The main conclusion that can be derived from this is that functions B and C are indeed non-convex, and thus require, as expected, a non-convex optimization method such as branch and bound to identify their global minimum. 8.2 Computing the Lower Bound In order to perform a search in a continuous search space, it is necessary to evaluate intervals of the control settings, instead of discrete values. In order to avoid confusion between control space intervals and time intervals, the 53
(a) Objective Function B (b) Objective Function C Figure 18: Histograms of objective function values of the local minima found for functions B and C with 100 random starting points, for a problem with PH = 7. time intervals i , with i = 1 ;:::;PH , will be referred to in this chapter as stages . The control actions v that have thus far been comprised of discrete values, are re-dened as a constrained control space box of NA dimensions: v = [ v 1 ;v 1 ] [ v 2 ;v 2 ] [ v NA ;v NA ] (45) Therefore, it is impossible to calculate the exact transition cost c of v . However, it is possible to determine the lowest possible value of c , given the intervals v of the actuators. This value can be used as the lower bound of the transition cost, denoted c L . Since the optimization problem makes use of an analytic model to calculate the state transition, as opposed to a black box model, it is well suited to apply the natural inclusion function to calculate its lower bounds [19]. This consists of exchanging the usual elementary operators by their interval extensions [19]: 8 > > > > < > > > > : [ a;b ]+[ c;d ] = [ a + c;b + d ] [ a;b ] [ c;d ] = [ a d;b c ] [ a;b ] [ c;d ] = [min f ac;ad;bc;bd g ; max f ac;ad;bc;bd g ] [ a;b ] [ c;d ] = [ a;b ] 1 d ; 1 c if 0 = 2 [ c;d ] The above given interval extensions were applied to each operation in the transition cost function c , in Eq.(40). It was hence also applied to the state transition function equations of the climate model. By implementing the extensions consistently throughout the equations, the natural inclusion function of the transition cost was obtained, denoted by F . The inclusion 54
function yields the maximum and minimum values of the transition cost c , given the actuator intervals v and state intervals s . F ( v ; s )=[ c ( v ; s ) ;c ( v ; s )] (46) The lowest value of the range was taken as the lower bound of the transition cost, c L : c L ( v ; s ) = min( F ( v ; s )) (47) The sum of c L up to i is a lower bound of the cost thus far: J L i = i X ` =1 c L ( v ` ; s ` ) 8 i = 1 ;:::;PH (48) It follows that the lower bound of the cost over a complete control sequence V with PH stages is: J L = J L PH = PH X i =1 c L ( v i ; s i ) (49) 8.3 Method 1: Multi-Stage Bisection This section describes the rst algorithm that was developed to search for a solution in a continuous search space. The aim of the developed algorithm is to stay close to the B&B search method developed thus far. To achieve this, the total control space of NA PH dimensions, was broken up into PH separate control space boxes of dimensions NA . In other words, each stage i has its own control space v i . At the start of the search, the box of each stage i is denoted v i 0 . The continuous set S 0 comprises the initial search boxes of all stages i = 1 ;:::;PH , and therefore encompasses the whole search space. S 0 = v 10 v 20 ::: v PH 0 Let the sub-set S t be a series of PH control boxes denoted v it . Therefore, S t is dened as: S t = v 1 t v 2 t ::: v PHt (50) 55
Figure 19: Graphical representation of individual control space for each stage i , narrowed down in each trajectory t with a bisection. To minimize the data storage, the algorithm performs a depth rst search. Therefore, instead of rening each control box v i in depth separately, the algorithm makes use of multiple trajectories, each going through all the control stages 1 to PH . Let a trajectory, indexed t , be a search through all stages i in the set S t 1 , yielding the subset S t . When the algorithm reaches the prediction horizon, it starts a new trajectory t +1, which searches the set S t by going through all the stages again and yields the set S t +1 . This process is performed repeatedly until the desired size of a subset is reached. The optimization problem outlined in Section 6 has three actuators (windows, fog system and heating), each with a control range of [0,100]. The set S 0 of control space boxes before the start of the search, at t = 0, is thus: 56
bound, disregarding the stage aspect of the problem, outlined in Section 4.1. This method was applied in order to have a benchmark algorithm for the other two continuous algorithms that were developed. The single-stage bisection searches the whole control space over all time stages at once. A problem with NA actuators and prediction horizon PH will have PH NA variables in its control sequence V . Each variable is denoted by v ai , with a = 1 ;:::;NA and i = 1 ;:::;PH . Recall Eq.(5): V = 0 B B B @ v 11 v 21 v PH 1 v 12 v 22 v PH 2 . . .. . ..... . . v 1 NA v 2 NA v PHNA 1 C C C A The search starts with a set S 0 enclosing the whole feasible area, stored in . The set is removed from and bisected into two subsets S 1 and S 2 by the longest edge. The algorithm calculates the lower bound J L for both subsets, according to Eq.(49). If J L > J U , the subset is discarded. In the case that J L < J U , the subset is stored in . Furthermore, the left sides of the subset intervals are taken as a feasible control sequence V f . If the objective cost J ( V f ) is lower than J U , then the upper bound is updated by J U = J ( V f ) and the optimal solution V = V f . Any subsets in with J L lower than the new upper bound are removed. After both subsets are either discarded or stored, the algorithm takes a new subset from and repeats the same steps as outlined for S 1 . As this process continues, the width of the sets diminishes. Sets that have a maximum width of are not bisected anymore, only evaluated and removed from . The branching and bounding process continues until is empty. The best solution thus far, V , hence becomes the best overall solution. The logical steps taken by the algorithm are outlined in Algorithm 7. 63
Algorithm 7 SingleStageBisection( S 0 ;; s 1 ) Require: S 0 ; and s 1 Determine upper bound J U Store S 1 in ; r := 1 while 6 = ? Remove a subset S from and split into 2 new subsets S r +1 and S r +2 Determine lower bounds J L r +1 and J L r +2 for p = r + 1 to r + 2 do if J L p < J U determine a feasible point V f 2 S p and J ( V f ) = J if J < J U J U := J remove all sets S k from with J L k > J U elseif Size ( S p ) > , Store S p in r := r + 2 endwhile return: V , J U 8.6 Comparison of Methods The alternative optimization methods outlined in Sections 8.3, 8.4, and 8.5 were programmed in Matlab 2015b in order to compare them. The accuracy of the methods was checked by comparing their solutions. Due to the determininstic nature of the optimization problem, all three methods should yield exactly the same control sequence V , down to the last decimal. Since this is the case, it is assumed that all three methods are correct. To inspect the quality of the methods further, their computational time was compared. The example problem of Section 6.2 was solved using all three alternative methods, at increasing prediction and control horizon lengths. Runs that took longer than 500 seconds were stopped and recorded at this value. The results are shown in Figure 22. 64
(a) Prediction Horizon (b) Control Horizon Figure 22: Computational time required by the three alternative solution methods Multi-stage Bisection, Multi-stage Multisection, and Single-stage Bisection, to solve the optimization problem described as example problem in Section 6.2 at increasing PH and CH values respectively. All runs were performed with = 1. All runs in (a) were performed with a CH : PH ratio of 1:3, and all runs in (b) had a prediction horizon of PH = 6. Figure 22a shows in increasing computational time with increasing prediction horizon PH length. This was expected, as the same pattern was found for the discrete algorithm in Figures 15 and 13d. Notice that for both the discrete and the continuous algorithm the computational time drastically increases at PH > 6. In the continuous control space, at runs with PH < 4 neither of the three methods shows an advantage over the others. For runs with 4 6 PH 6 6, the single-stage bisection has the lowest computational time, followed by the multi-stage bisection, while the multi-stage octosection is the slowest. The same order of computational speed occurs in Figure 22b. The graph also shows that the maximum control horizon that the algorithms can handle within the 500s time span, given PH = 6, is CH = 2. The computational time appears to increase exponentially with CH , since the R 2 values of the exponential trendlines are all close to 1. The reason for this pattern is that, since the branching only occurs before the control horizon, for a longer control horizon CH the complexity of the problem increases. Overall, all three methods can be concluded to be able to handle only relatively low complexity problems. For the multi-stage multisection, the reason for the high computational time is the large number of node evaluations that must be performed. Since the control space at each stage i is split 65
in eight, the search tree grows very wide and only small sections of the tree can be pruned. The multi-stage bisection splits the control space in two. Therefore, it can bound larger tree sections and achieve a lower computational time than the octosection. However, since the search space is merely bisected, three times as many splittings are needed to achieve the same interval width . This method therefore requires a much deeper search. In the single-stage bisection method, the high computational time occurs because the state and control intervals of the inclusion function widen at every stage i , hence making the lower bounds less accurate. Consequently, for a larger prediction horizon PH , the range of the inclusion function becomes too wide and the lower bounds are not high enough to prune large sub-trees. Regarding the applicability of the methods, it is necessary to recall from Section 5.1 that time was discretized into intervals of 15 s, each of which is one stage. In MPC, one optimization problem must be solved at every stage. Therefore, the computational time of one optimization problem may not be higher than 15 s. However, the computational time is dependent on the implementation medium and platform. As a consequence, conclusions can only be made about the computational time of the methods with respect to each other. The absolute values of the computational time could be reduced by running them on other programs, so it is not possible to conclude whether the algorithms would be applicable in reality. In both graphs, the single-stage bisection is faster than the multi-stage searches. Consequently, it can be concluded that the branch and bound application from literature does not have any advantages over a regular branch and bound search when applied in continuous control space. 66
9 Discussion and Conclusion 9.1 Discussion This thesis investigated the branch and bound algorithm applied to MPC with variants for a discrete and continuous control space. The results of the experiments were already discussed in previous chapters; therefore, this chapter focuses on a discussion of the assumptions that were made throughout the report. To evaluate results of the experiments, the computational time of the algorithms was compared. However, it must be kept in mind that computational time can also be inuenced by possible programming ineciencies, which are caused by the programmer and are hence unrelated to the algorithm itself. On the other hand, algorithms can also be compared using the number of node evaluations. Yet if the algorithms that are being compared do not perform the same computation in a node, the number of node evaluations is also not a valid method of comparison. Therefore, this leads to the discussion point of the validity of comparisons. In order to make a precise comparison, it would be necessary to revise the code of the models to ensure that the results are not inuenced by the programmer. Another important issue is the validity of the objective functions that were implemented. Although the functions were taken from literature, they are a simplistic approach to describing the objectives of the user. In reality, climate control is a multi-objective problem, and should hence be tackled as such. This puts the results of the report in question, as it raises the doubt whether the branch and bound algorithm would still be applicable and useful in the case of a multi-objective optimization problem. Finally, the nature of the optimization problem leads to the notion that it might also be possible to solve it using dynamic programming. This technique also uses a transition function to optimize variables in time. However, dynamic programming searches from the last to the rst time stage, as opposed to the branch and bound method, which searches the stages in their chronological order. To further evaluate the quality of the branch and bound method, it would be of interest to compare it to dynamic programming to determine which of the two is more suitable. 9.2 Conclusion The main purpose of this study was to determine the applicability of branch and bound to climate control in greenhouses. Initially, it was found that climate control is usually performed using Proportional Integral Derivative (PID) controllers, which correct for past disturbances of the climate based on simplied transfer function models [28]. An alternative control method is Model Predicitive Control (MPC), which controls the actuators based on a prediction of the future climate state. MPC discretizes time into intervals 67
and repeatedly optimizes a control sequence V over a prediction horizon of PH time intervals. Literature shows that branch and bound can be used as an optimization method for the control sequence. The B&B optimization of control sequence V was examined in the available literature. It was found that the algorithm calculates the transition cost of a discrete control action v at every time interval of the prediction horizon, and bounds the search if the sum of the transition costs up to that interval is higher than an upper bound. To test the applicability of the search method, it was extended to a model with climate simulation. The literature was reviewed regarding objective functions for climate control. A comparison was made of three alternative functions by implementing them in the previously mentioned model. It was found that the most suitable objective function is one that does not aim to follow a specic setpoint, but instead keeps the climate state variables within a certain range. By implementing a penalty and an oset for solutions violating the allowed range, the function reduces the computational time while adhering to the desired climatic conditions in the greenhouse. From the implementation of the programmed model it was concluded that the B&B search could prove a useful optimization method. Its main advantages are that it always nds the global optimum, it deals with constraints implicitly, and it requires no initial guess of the solution. The computational time of the method was found to be mostly dependent on the prediction horizon and the outside temperature an input parameter of the climate model. At increasing prediction horizon PH lengths the computational time increases exponentially, while increasing outside temperatures within the permitted range leads to a linear computational time increase. Computational time for outside temperatures below and above the permitted range were found to be constant, yet higher for temperatures above the range. Overall, the B&B search was deemed applicable, given the specied parameter values, at prediction horizon of lengths lower than 6. The main disadvantage of the B&B optimization is that it requires a discretization of the control space. Consequently, it requires a trade-o between discretization and performance [20]. A possible computational time reduction was investigated through a modication of the algorithm. The modied algorithm evaluates all intervals after the control horizon CH in one single node by simulating all the states until the end in one step. This was possible because after the control horizon all control actions stay constant. The modication was tested in a programmed model, and found to largely reduce computational time of optimization problems with high computational eort. It was therefore concluded that the modication could be a promising addition to the algorithm. Finally, the possible use of a continuous control space was examined. The control space was analyzed and concluded to have multiple minima, hence justifying the implementation of a B&B search. The natural extension 68
function of the objective function was computed using interval arithmetics, to calculate lower bounds of a subset. Three alternative search methods were developed. Multi-stage bisection adheres to the B&B implementation used thus far, by maintaining the principle of a cumulative objective value throughout the time stages. Each time stage has its own control space, which is bisected repeatedly until a certain accuracy degree . Multi-stage multisection applies the same principle as multi-stage bisection, but bisects all edges of a control space, rather than just the longest edge. Finally, the single-stage bisection performs a regular B&B search, by bisecting the whole control space at once. After implementing the methods in a programmed model, all three were concluded to have an exponentially increasing computational time at increasing prediction as well as control horizon. The single-stage bisection was deemed the most suitable method, since it could handle the same PH and CH values with lower computational times. This leads to the conclusion the B&B algorithm from literature does not show any advantages when applied in continuous control space. Overall, the main conclusions of the investigation is that B&B as an optimization method for MPC in a greenhouse is a promising search method if applied in discrete control space, its computational time can be reduced by simulating the last nodes of the control space in one step, and it does not show advantages in continuous control space. 9.3 Further Research The foremost topic requiring further research in this topic is the selection rule that was applied when branching. In the models programmed for this report, the B&B algorithms were programmed to choose branches from left to right. However, it could prove benecial to use a more elaborate selection method, for example, based on the transition cost, to choose the most promising branch. This could help to reach a lower upper bound faster, and hence prune many more sections of the tree. Another topic for further research is the application of the modication for node reduction to the continuous problem. Since the modication led to a reduction in computational time in the discrete control space, the node reduction could have the same eect on the computational time of the continuous search algorithms. In the continuous search space it would also be of interest to investigate the bisection technique that is applied to the control spaces. In this report, the bisection was always performed on the largest edge of a control set, as this is a common technique. However, literature has shown that largest edge bisection does not always yield the smallest search tree [31]. It would thus be of interest which bisection technique would be optimal for this particular search algorithm, and whether a dierent technique could improve its applicability in continuous search space. Finally, in order to determine the applicability of the B&B algorithm 69
in MPC with more certainty, it would be necessary to compare its performance to that of other optimization techniques. In particular dynamic programming would be of interest to compare with, as this method also optimizes problems over multiple stages or intervals. Finally, the B&B algorithm should be put in a broader context, to investigate its applicability to control problems in other elds. 70
References [1] J. M. Aaslyng, J. B. Lund, N. Ehler, and E. Rosenqvist. Intelligrow: a greenhouse component-based climate control system. Environmental Modelling & Software , 18(7):657{666, 2003. [2] J. Baldwin and B. Pilsworth. Dynamic programming for fuzzy systems with fuzzy environment. Journal of mathematical analysis and applications , 85(1):1{23, 1982. [3] W. Baudoin, R. Nono-Womdim, N. Lutaladio, A. Hodder, N. Castilla, C. Leonardi, S. De Pascale, and M. Qaryouti. Good agricultural practices for greenhouse vegetable crops: principles for mediterranean climate areas . Number 217. Food and Agriculture Organization of the United Nations, 2013. [4] M. Berenguel, F. Rodrguez, F. Acien, and J. Garca. Model predictive control of ph in tubular photobioreactors. Journal of Process Control , 14(4):377{387, 2004. [5] X. Blasco, M. Martnez, J. Herrero, C. Ramos, and J. Sanchis. Modelbased predictive control of greenhouse climate for reducing energy and water consumption. Computers and Electronics in Agriculture , 55(1):49{70, 2007. [6] N. Castilla and F. Nuez. Manejo del cultivo intensivo con suelo. El cultivo de tomate. Madrid: Mundi-prensa , pages 212{217, 1995. [7] P. Ferreira and A. Ruano. Discrete model based greenhouse environmental control using the branch & bound algorithm. In Proceedings of the 17th IFAC World Congress, Seoul, Korea , volume 611, 2008. [8] P. M. Ferreira and A. E. Ruano. Application of computational intelligence methods to greenhouse environmental modelling. In Neural Networks, 2008. IJCNN 2008.(IEEE World Congress on Computational Intelligence). IEEE International Joint Conference on , pages 3582{3589. IEEE, 2008. [9] M. C. Garca-Martnez, S. Balasch, F. Alcon, and M. FernandezZamudio. Characterization of technological levels in mediterranean horticultural greenhouses. Spanish Journal of Agricultural Research , 8(3):509{525, 2010. [10] P. E. Gill, W. Murray, and M. H. Wright. Practical optimization . Academic press, 1981. 71
[11] F. He, C. Ma, J. Zhang, and Y. Chen. Greenhouse air temperature and humidity prediction based on improved bp neural network and genetic algorithm. In Advances in Neural Networks{ISNN 2007 , pages 973{980. Springer, 2007. [12] E. M. T. Hendrix and B. G.-. Toth. Introduction to Nonlinear and Global Optimization . Springer, New York, 2010. [13] J. Herrero, X. Blasco, M. Martnez, C. Ramos, and J. Sanchis. Nonlinear robust identication of a greenhouse model using multi-objective evolutionary algorithms. Biosystems Engineering , 98(3):335{346, 2007. [14] E. Heuvelink. Inuence of day and night temperature on the growth of young tomato plants. Scientia Horticulturae , 38(1-2):11{22, 1989. [15] E. Horowitz and S. Sahni. Fundamentals of computer algorithms . Computer Science Press, 1978. [16] R. Johansson. System modelling and identication. Automatica , 3(33):477{478, 1997. [17] H.-P. Liebig. Temperature integration by kohlrabi growth. In Symposium on High Technology in Protected Cultivation 230 , pages 371{380, 1988. [18] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. Scokaert. Constrained model predictive control: Stability and optimality. Automatica , 36(6):789{814, 2000. [19] I. Mazhoud, K. Hadj-Hamou, J. Bigeon, and G. Remy. Interval-based global optimization in engineering using model reformulation and constraint propagation. Engineering Applications of Articial Intelligence , 25(2):404{417, 2012. [20] L. F. Mendonca, J. Sousa, and J. S. da Costa. Optimization problems in multivariable fuzzy predictive control. International Journal of Approximate Reasoning , 36(3):199{221, 2004. [21] L. Miranda, I. Schuch, D. Dannehl, T. Rocksch, R. Salazar, and U. Schmidt. Using articial neural networks to predict the climate in a greenhouse: First simulation results on a semi-closed system. In II International Symposium on Horticulture in Europe 1099 , pages 137{144, 2012. [22] L. Mitten. Branch-and-bound methods: General formulation and properties. Operations Research , 18(1):24{34, 1970. [23] J. A. Nelder and R. Mead. A simplex method for function minimization. The computer journal , 7(4):308{313, 1965. 72
Crop evapotranspiration: G 3 ( s 1 ;s 2 ) = 2 ( F 8 ( s 1 ) 26 + 2 14 24 5 F 9 ( s 1 ;s 2 ) 8 ) ( F 8 ( s 1 ) + 23 (1 + ( 8 =F 12 ( s 1 ;s 2 )))) F 10 ( s 1 )(A.11) F 8 ( s 1 ) = f 1 ( s 1 + 0 : 5) f 1 ( s 1 0 : 5) (A.12) 26 = (1 e 12 14 ) 25 d 3 (A.13) F 9 ( s 1 ;s 2 ) = F 3 ( s 1 ) 1 s 2 100 (A.14) F 10 ( s 1 ) = (3 : 1468 0 : 002365( s 1 + 273)) 10 6 (A.15) F 11 ( s 1 ;s 2 ) = ( 0 : 39 0 : 029+ F 9 ( s 1 ;s 2 ) ; F 9 ( s 1 ;s 2 ) 0 : 361 1 ; F 9 ( s 1 ;s 2 ) < 0 : 361 (A.16) F 12 ( s 1 ;s 2 ) = 10 + ( 9 10 ) 1 exp 25 d 3 160 F 11 ( s 1 ;s 2 ) (A.17) (A.18) Solar energy supplied to air volume: Q 2 ( d 3 ) = 2 25 d 3 (A.19) Energy exchange by conduction and convection phenomena: Q 3 ( s 1 ;d 1 ) = 2 15 ( s 1 d 1 ) (A.20) Energy loss due to crop evapotranspiration: Q 4 ( s 1 ;s 2 ) = F 10 ( s 1 ) G 3 ( s 1 ;s 2 ) (A.21) Energy exchange due to window ventilation: Q 5 ( v 1 ;s 1 ;d 1 ) = 24 5 F 2 ( v 1 )( s 1 d 1 ) (A.22) Energy loss by nebulization: Q 6 ( v 2 ;s 1 ) = F 10 ( s 1 ) G 1 ( v 2 ) (A.23) Energy exchange between thermal mass and inside air: Q 7 ( s 1 ;s 3 ) = 2 11 ( s 3 s 1 ) (A.24) Energy stored by the thermal mass during the day: Q 8 ( d 3 ) = 21 Q 2 ( d 3 ) (A.25) Energy loss through ground: Q 9 ( s 3 ) = 2 13 s 3 17 20 (A.26) 79
Functions Saturation pressure: f 1 ( Temp )=0 : 61(1 + 1 : 414 sin (5 : 82 e 3 s 1 emp )) 8 : 827 (A.27) Relative to absolute humidity: f 2 ( Temp;Rel:H ) = Rel:H 0 : 611 F 3 ( s 1 ) 100 16 (A.28) Absolute to relative humidity: f 3 ( Temp;Abs:H ) = ( 100 ; f 4 ( Temp;Abs:H ) > 100 f 4 ( Temp;Abs:H ) ; f 4 ( Temp;Abs:H ) 100 (A.29) f 4 ( Temp;Abs:H ) = ( Abs:H )100 16 0 : 611 f 1 ( Temp )(A.30) 80