Robust solutions of bi-blend recipe optimization with quadratic constraints

Robust solutions of bi-blend recipe optimization with quadratic constraints Author: Freek Niewerth MSc Research report Universidad de M´alaga in co-operation with Wageningen University and Universidad de Almer´ıa Dpt. Computer Architecture, Operations Research and Logistics and Computer Architecture and Electronics Supervisors: E.M.T. Hendrix and L.G. Casado February 26, 2009 2 Abstract Production companies use raw materials to compose end-products. They often make different products with the same raw materials. In this research, the focus lies on the production of two end-products consisting of (partly) the same raw materials as cheap as possible. Each of the products has its own demand and quality requirements consisting of quadratic constraints. The minimization of the costs, given the quadratic constraints is a global optimization problem, which can be difficult because of possible local optima. Therefore, the multi modal character of the (bi-) blend problem is investigated. Standard optimization packages (solvers) in Matlab and GAMS were tested on their ability to solve the problem. In total 20 test cases were generated and taken from literature to test solvers on their effectiveness and efficiency to solve the problem. The research also gives insight in adjusting the quadratic constraints of the problem in order to make a robust problem formulation of the bi-blend problem. Contents 1 Introduction 4 1.1 Generalgoal............................... 4 1.2 Globaloptimization........................... 4 1.3 Problemstatement ........................... 6 1.4 Vision .................................. 6 1.5 Researchquestions ........................... 6 1.6 Structure of this report . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 Blending with quadratic constraints 7 2.1 Unit simplex and Cost function . . . . . . . . . . . . . . . . . . . . 7 2.2 Linearconstraints............................ 8 2.3 Quadratic constraints . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4 Formulation of blending problem . . . . . . . . . . . . . . . . . . . 10 2.5 Global and local optima . . . . . . . . . . . . . . . . . . . . . . . . 11 2.5.1 Feasible area with two compartments . . . . . . . . . . . . . 12 2.5.2 Non-convex feasible compartment . . . . . . . . . . . . . . . 12 2.6 Testcases ................................ 13 3 Bi-Blending 15 3.1 Introduction............................... 15 3.2 Capacity constraint . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.3 Costfunction .............................. 16 3.4 Formulation of the Bi-blend problem . . . . . . . . . . . . . . . . . 16 3.5 A 2-dimensional example . . . . . . . . . . . . . . . . . . . . . . . . 16 3.6 Testcases ................................ 19 4 Solvers 21 4.1 Localsolvers............................... 21 4.1.1 Multistart............................ 22 4.2 Globalsolvers.............................. 23 4.3 Performance indicators . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.3.1 Effectiveness........................... 23 4.3.2 Efficiency ............................ 23 4.4 Blendingcases.............................. 24 4.4.1 Experimental specifications . . . . . . . . . . . . . . . . . . 24 4.4.2 Localsolutions ......................... 24 2 Contents 4.4.3 Performance........................... 25 4.4.4 Conclusion............................ 27 4.5 Bi-blendcases.............................. 28 4.5.1 Localsolutions ......................... 28 4.5.2 Performance........................... 29 4.5.3 Conclusion............................ 29 5 Robustness 31 5.1 Introduction............................... 31 5.2 -Robustness with respect to (bi-)blending . . . . . . . . . . . . . . 32 5.3 Some (heuristic) approaches . . . . . . . . . . . . . . . . . . . . . . 32 5.3.1 Reformulation of quadratic constraints . . . . . . . . . . . . 33 5.3.2 -robustness on X........................ 37 5.4 Conclusion................................ 38 6 Conclusion 40 A Testcases 42 A.1 Quadratic constraints . . . . . . . . . . . . . . . . . . . . . . . . . . 45 B GAMS code 55 B.1 Rumcoke................................. 55 B.2 Bi-blendcase1 ............................. 57 C Matlab code 59 C.1 Rumcoke................................. 59 C.1.1 f(x) ............................... 59 C.1.2 g(x) ............................... 59 C.1.3 Fmincon single start . . . . . . . . . . . . . . . . . . . . . . 59 C.1.4 Fmincon-multi ......................... 60 Bibliography 61 3 Chapter 1 Introduction 1.1 General goal This master thesis is written for Wageningen University and the University of Almeria. In order to write a good report the University of Wageningen has set some predefined standards. These standards are described in the Msc Thesis protocol of the WUR. The general goal of a thesis written by a WUR-student is the following: “The overall goal of the thesis is the development of research skills and the ability to analyze and present research results in a systematic and clear way. The thesis is the culmination of the MSc study program in which the student will have to show that he/she is able to design and conduct social science research at an academic level and is able to theoretically reflect on a particular field of research relevant to the MSc program at hand.” [1] 1.2 Global optimization This thesis is about a specific global optimization problem. The objective of a global optimization method is to find the best possible solution (global optimum) of a problem by modifying decision variables in order to minimize (or maximize) an objective function, while there may also exist some local optima. Sometimes the decision variables have to fulfil some requirements. This can be compared with the search of the highest point above sea level in the world. When we want to maximize the height (objective) of mountains in the world we can vary the location by adjusting the longitude and the latitude. These are the decision variables for this problem. A local optimum in the region of the Alps will be Mont Blanc, in the region of Africa this will be the Kilimanjaro and in the region of Asia the Mount Everest. The global optimum however is the Mount Everest since this is the highest mountain in the world. In Table 1.1 the main characteristics of optima are described. We can see that the local optimum for Asia is also the global optimum for the world. This gives us the following property: Global optimum ⊂Local optima. 4 Chapter 1: Introduction Global optimization problems can be difficult because of the existence of local optima. When one would only search in the Alps and conclude Mont Blanc is highest in that region, one could easily think this is the highest mountain of the world since in its direct neighbourhood there is no mountain which is higher than Mont Blanc. Table 1.1: Optima Type Corresponding mountain Global optimum Mount Everest Local optimum Mont Blanc, Kilimanjaro, Mount Everest This thesis is about a global optimization problem in which minimizing a function given some decision variables is the objective. We will use the global optimization problem formulation from [5] to describe optimization problem in a mathematical way as is shown in (1.1). min{f(x)}, x ∈X⊂Rn,(1.1) in which f(x) is a real valued continuous function and xvaries in a continuous way in a set Xwith a dimension n. Formula (1.1) describes the minimization of a function f. We are trying to find the global minimum, where there also may exist local minima. In Figure 1.1, function f(x) is plotted of which we want to know the minimum value on the interval 0 ≤x≤20. In the figure we can see four local optima that are minimal in their environment. The lowest local optimum is the global optimum and the value we want to know. -3 -2 -1 0 1 2 3 4 5 6 7 0 5 10 15 20 Global optimum Local optima Local optimum Figure 1.1: Four local optima for f(x) = 2 cos(x)−x 4+ 4, 0 ≤x≤20. In this thesis, the global optimization problem of mixing raw ingredients in order to produce two different products will be discussed; the so-called bi-blending problem. 5 Chapter 1: Introduction 1.3 Problem statement Companies use raw materials to make products. In fact, they often make different products with the same raw materials. In this research we focus on a problem with two products consisting of (partly) the same raw materials. Each of the products has its own demand and quality requirements consisting of quadratic constraints. In order to optimize this so called bi-blending problem we need to find robust solutions. To decide whether a solution is robust the following definition of a robust decision is used: “A robust decision is the best possible choice, one found by eliminating all the uncertainty possible within available resources, and then choosing, with known and acceptable levels of satisfaction and risk” [11]. This means for example that a found optimal solution must stay feasible despite some known variation in dosage of a raw material. 1.4 Vision The main goal of the project is to investigate methods to optimize and find robust solutions of bi-blending problems with quadratic constraints. We need to develop a way to generate (robust) solutions numerically. This can either be done with existing software and/or by designing new algorithms. The vision of this thesis is to find a standard procedure to optimize biblending problems in a robust way. 1.5 Research questions Corresponding to the problem statement and vision the main question of this thesis is: •Which methods can be used to optimize and find robust solutions of biblending problems with quadratic constraints? This question is subdivided in four sub questions: 1. How can the optimization problem be described? 2. What is the multimodal character of the problem? 3. What is the potential of existing standard software? 4. How can the robustness of solutions be determined? 1.6 Structure of this report In Chapter 2, the blending problem is discussed, in Chapter 3 the bi-blending problem. Chapter 4 gives results of standard optimization software on (bi-) blending problems and in Chapter 5 the robustness problem is further investigated. 6 Chapter 2 Blending with quadratic constraints In the food industry, raw materials are put together and processed to produce products. The easiest way of processing is done by putting different raw materials together to mix them. These mixing processes do not only occur in the food sector but also in other types of industries. 2.1 Unit simplex and Cost function The blending of raw materials can be described in a mathematical way. A vector xwith a number of elements equal to the amount of raw materials is used to indicate the fraction (xi) of every raw material (i) used in the mix. All fractions in the mix have to add up to 1. A recipe (of a mix) is mathematically defined by the unit simplex: S={x∈Rn| n X i=1 xi= 1; xi≥0}.(2.1) where ndenotes the number of raw materials. In Figure 2.1b, the unit simplex is given in R3(three raw materials). In this simplex we can find 7 possible sets depending on the composure and the number of raw materials involved. Namely, a recipe consisting of only: x1,x2or x3, denoted by the points; a recipe consisting of: x1and x2,x1and x3or x2and x3, denoted by the lines; and a recipe consisting of x1,x2and x3, denoted by the triangular area. 7 Chapter 2: Blending with quadratic constraints using the 3 dimensional unit simplex from Figure 2.2. For the 3-dimensional cases we have the ’Rum-coke’ problem and 10 other test cases. Test case 1 to 6 are ‘regular’ blending problems with a feasible area which consists of one compartment or where no feasible solution exists. Test cases 7 to 10 are extreme cases. Test case 7 is the non-convex problem from Figure 2.7; it may have multiple local optima when an appropriate cost vector cis chosen. Test case 8 to 10 each have two compartments for their feasible area and thus certainly contain local optima, like Example 2.3. Although the plot in Figure 2.5 seems straightforward in finding a feasible area, solving a blending problem is not easy. Making a plot can be time consuming and for dimensions higher than 3 a plot cannot be made. Therefore, algorithms are used for solving quadratic blending problems. In some software, like Excel, Matlab and GAMS, solvers exist which can help us solving the problems. In Chapter 4 several of these solvers will be tested and compared with each other. In the next chapter the concept of bi-blending is explained. In many industries, several products are made simultaneously with use of the same raw materials, which creates a new problem when certain raw materials become scarce. 14 Chapter 3 Bi-Blending 3.1 Introduction The focus of this research is on bi-blending. The concept of bi-blending is taken from multi-blending where multiple products are made from the same (scarce) raw materials. For this research, the focus is on two products; this is called the bi-blend problem. In industry manufacturers sometimes face the problem of scarcity of raw materials for the products they want to make. As with the blending problem, this can be solved by adding a bigger dosage of another ingredient. For the description of the two recipes, we use the variables xand yfor product 1 and product 2 respectively. The scarcity of raw materials is described by a capacity constraint. 3.2 Capacity constraint The availability of raw materials is described by the capacity constraint. This constraint is given by (3.1). D1xi+D2yi≤Bii= 1, . . . , n, (3.1) where D1and D2represent the amount of respectively product 1 and product 2 to be made. The amount of available raw material iis given by Bi. Example 3.1 A farmer wants to grow 10 pigs and 3 cows. A pig needs fodder consisting of 80% corn (i= 1) and 20% water (i= 2) and a cow respectively 70% and 30%. However, one cow uses twice the amount of raw materials. The farmer has a silo with a capacity of 15 units of corn. Water is unlimited. The farmer wants to know if his silo is big enough to grow the pigs and the cows. When we enter the values for the variables and data (x, y, D and B)we can easily decide if the silo is big enough: x=0.8 0.2, y =0.7 0.3, D1= 10, D2= 3 ∗2=6, B =15 ∞ . 15 Chapter 3: Bi-Blending Substitution in (3.1) gives: 10 ∗0.8 0.2+ 6 ∗0.7 0.3≤15 ∞ . The sum of the used corn (i= 1) is: 10 ∗0.8+6∗0.7 = 12.2≤15. The silo is big enough so there are no problems for the farmer. 3.3 Cost function With the capacity constraint added to the problem, the cost function changes since we are now taking two products into account with their individual demands. The cost function of the bi-blending problem can be written as: f(x, y) = n X i=1 ci(D1xi+D2yi).(3.2) We can add the capacity constraint and the new cost function to quadratic constraints for each individual product to formulate the bi-blend problem. 3.4 Formulation of the Bi-blend problem The Bi-Blending problem is described by (3.3). min {f(x, y) = Pn i=1 ci(D1xi+D2yi)}(Cost (3.2)) s.t. x, y ∈S(Blending (2.1)) x∈Q1, y ∈Q2(Feasibility (2.3)) D1xi+D2yi≤Bii= 1, . . . , m (Capacity constraint (3.1)). (3.3) Note that the ingredients are denoted by i. This means that xiand yi(for the same i) are fractions of the same ingredient; only the amount of the ingredient varies. We give an easy example of a bi-blend problem with only bounds on the use of raw materials for product 1 and product 2 and a capacity constraint. 3.5 A 2-dimensional example We have an instance with bounds on the mixtures for product 1 and 2. The parameters for this problem are: c=1 2, D1= 2, D2= 4, x ≤0.8 0.4, y ≤0.4 0.8, B =5 5. The problem is sketched in Figure 3.1. 16 Chapter 3: Bi-Blending 0 0.2 0.4 0.6 0.8 1 1.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 x1 x2 X* (a) x* 0 0.2 0.4 0.6 0.8 1 1.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 y1 y2 Y* (b) y* Figure 3.1: Bi-blend problem 1 Table 3.1: Data and solution of Bi-blend instance 1 i xiD1yiD2T otali Bi 1 0.8 2 0.4 4 3.2 5 2 0.2 2 0.6 4 2.8 5 f* = 3.2 * 1 + 2.8 * 2 = 8.8 We can see that the capacity constraint is not binding in Figure 3.1 and from Table 3.1. Of raw material i=1, 2*0.8 + 4*0.4 = 3.2 is used and of raw material i=2, 2*0.2 + 4*0.6 = 2.8 is used while the capacity is BT= (5,5). In both product 1 and 2 the maximum amount of i= 1 is used since this is the cheapest ingredient. This means that for product 1 and 2 the cheapest mixtures could be made. This means that when the capacity constraint is not binding we can write the bi-blend problem as two separate blending problems as in (3.4). min x,y∈Q∩S{f(x, y) = n X i=1 ci(D1xi+D2yi)}= min x∈Q1∩S{D1cTx}+ min y∈Q2∩S{D2cTy}(3.4) When the capacity constraint is binding, the problem is more interesting. Lets set BT= (3,5) (so that B1≤3.2 since i= 1 is the cheapest ingredient) and observe what happens. Figures (and corresponding Tables) 3.2 and 3.3 give local optima that are extreme points for this example. 17 Chapter 3: Bi-Blending 0 0.2 0.4 0.6 0.8 1 1.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 x1 x2 X* (a) x* 0 0.2 0.4 0.6 0.8 1 1.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 y1 y2 Y* (b) y* Figure 3.2: Extreme point 1 for Bi-blend problem 1 Table 3.2: Extreme point 1 for Bi-blend instance 1 i xiD1yiD2T otali Bi 1 0.7 2 0.4 4 3 3 2 0.3 2 0.6 4 3 5 f*=3*1+3*2=9 0 0.2 0.4 0.6 0.8 1 1.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 x1 x2 X* (a) x* 0 0.2 0.4 0.6 0.8 1 1.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 y1 y2 Y* (b) y* Figure 3.3: Extreme point 2 for Bi-blend problem 1 Table 3.3: Extreme point 2 for Bi-blend instance 1 i xiD1yiD2T otali Bi 1 0.8 2 0.35 4 3 3 2 0.2 2 0.65 4 3 5 18 Chapter 3: Bi-Blending f*=3*1+3*2=9 From Figures 3.2 and 3.3 we can conclude that there are infinitely many local optima. This is caused by the binding capacity constraint on the cheapest ingredient (i= 1). The sum of x1and y1will always be 3 (and accordingly the sum of x2and y2will always be 3), but these ingredients can be distributed among product 1 and 2 according to their individual linear constraints in any way we want. The above mentioned local optima are boundaries of possible recipes. The set of global optimal solutions is a line piece. More alternative solutions are given in Table 3.4. When we have a multidimensional problem, we can get a set of global solutions which covers a multidimensional area when all most cheapest ingredients are binding on the capacity constraint and there exists a feasible solution. Table 3.4: Alternative solutions for Bi-blend instance 1 i xiD1yiD2T otali Bi 1 0.8 2 0.35 4 3 3 2 0.2 2 0.65 4 3 5 1 0.775 2 0.3625 4 3 3 2 0.225 2 0.6375 4 3 5 1 0.75 2 0.375 4 3 3 2 0.25 2 0.625 4 3 5 1 0.725 2 0.3875 4 3 3 2 0.275 2 0.6125 4 3 5 1 0.7 2 0.4 4 3 3 2 0.3 2 0.6 4 3 5 f*=3*1+3*2=9 3.6 Test cases We construct bi-blend cases using quadratic constraints from Appendix A.1 to make restrictions on product 1 and product 2. We add a cost function and a capacity constraint. The demand Din all these cases will be DT= (1,1). We made five 3-dimensional bi-blend test cases which are described in Table 3.5. In every row of the table a test case is given with their corresponding quadratic constraints on xand yand the used cost function and capacity constraint. 19 Chapter 3: Bi-Blending Table 3.5: Bi-blend test cases test case Q.constraints on xQ.constraints on yCost function Capacity vector BT 1 1, 2 22, 24 (1.1, 1.7, 2.0) (2, 2, 2) 2 1, 2 22, 24 (1.1, 1.7, 2.0) (0.5, 2, 2) 3 1, 2 22, 24 (1.1, 1.7, 2.0) (0.5, 0.75, 2) 4 1, 10 11, 15 (3.0, 2.0, 1.0) (2, 2, 2) 5 1, 10 11, 15 (3.0, 2.0, 1.0) (2, 2, 0.7) 6 1, 10 11, 15 (3.0, 2.0, 1.0) (2, 0.7, 0.7) 7 22, 24 2, 6 (3.0, 2.0, 1.0) (2, 2, 2) 8 22, 24 2, 6 (3.0, 2.0, 1.0) (2, 2, 0.5) 9 22, 24 2, 6 (3.0, 2.0, 1.0) (2, 0.9, 0.5) 10 22, 24 2, 6 (3.0, 2.0, 1.0) (0.42, 2, 0.73) In the next chapter, optimization algorithms are tested on the test cases we made on (bi-) blending. 20 Chapter 4 Solvers There is no use of programming a new algorithm when there already exists a good optimization algorithm (solver) that can get the same results. In this chapter the potential of some existing solvers is examined. This chapter distinguishes two types of solvers; solvers that find local optima and solvers that find the global optimum. Instances of the (bi-) blending problem will be examined with different non-linear programming (NLP) solvers in GAMS and Matlab. In Table 4.1 the used solvers are given with their corresponding access language. Also the type of solver is indicated, which will be further explained in this chapter. There are more standard software packages which contain NLP solvers we did not consider. Table 4.1: Solvers Solver Access language Type of solver BARON GAMS Global LGO GAMS Global OQNLP GAMS Multi start Local MINOS GAMS Local SNOPT GAMS Local CONOPT GAMS Local Fmincon Matlab Local Fmincon-multi Matlab Multi start Local 4.1 Local solvers Many standard solvers search for local optima. This means that given a starting point the solver tries to converge to a local optimum using a local search method. A local search may return a local optimum where there may also exist better optima. These solvers are tested to observe how difficult (bi-) blending problems are, and whether they might be solved by a local solver. 21 Chapter 4: Solvers 4.1.1 Multi start In the case we want to find the global optimum of an instance while using a local solver, several local searches can be applied with new starting points for each search. In [9] a multi start method is used by placing a grid over the unit simplex to perform a local search from each point in the grid. Another way of doing many local searches is to use randomly chosen starting points. In this research this stochastic multi start method is applied in GAMS and Matlab. We can try to find all local optima using the multi start method and pick the best; however this does not guarantee to find the global optimum. “If after some calculation time no solution of the inequality problem has been found, it is not certain whether there exists one [5].” If on the other hand many local searches are done, the chance of not finding an existing global optimum becomes very small. With use of a loop a multi start can be done while generating the starting points randomly for each new iteration (a local search) of the loop. In Matlab this is easy using a for-loop with randomly chosen starting points on the unit simplex (2.1). We can generate uniformly distributed starting points on the unit simplex using Algorithm (1) from [10]. For Fmincon-multi, 200 random starting points were chosen per instance. Figure 4.1 shows the 3-dimensional unit simplex with the notation used in Figure 2.4 with 200 randomly generated starting points. Algorithm 1 : Generating uniform random points on the unit simplex Funct S(n) 1. for (i=1:n) ai∼Ne(1) 2. for (i=1:n) Si=ai Pn j=1 aj 3. return S 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X1 X2 X3 Figure 4.1: 200 uniformly distributed starting points on the unit simplex. 22 Chapter 4: Solvers In GAMS it is not straightforward to implement a for-loop, That is why no manual loop is implemented in GAMS. However, in GAMS exists a stochastic multi start solver named OQNLP in which the local solvers of GAMS chosen for this research can be implemented. A disadvantage in using this multi start algorithm is that it generates random starting points over a box, where we would like to have random starting points which are on the unit simplex. We can set bounds on the starting points so it will be near the unit simplex with 0 ≤x≤1 to give a solution to this problem. 4.2 Global solvers Global solvers are solvers which use information from the whole search region to find the global optimum. With this global information the global character of a found optimum can be guaranteed. In this research GAMS/BARON and GAMS/LGO are solvers which claim to be global optimizers using different methods. This means that they should return the global optimum if a problem has one. 4.3 Performance indicators 4.3.1 Effectiveness The performance of a solver can be measured with different indicators. The most important indicator for solvers is whether the solver is able to solve a given blending problem; this is called the effectiveness. The way solving ability is hard to measure. For instance a (bi-) blending problem can be feasible or infeasible. Before we can say a problem has a feasible solution or not, we have to find one or we have to prove there is none. If we do find a feasible optimum with a given solver we do not know whether it is global unless the solver can prove it. Only global solvers can make this guaranty, for example by exploring lower bounds over the whole search region. Multi-start local solvers can be very good as well in finding a global optimum, but where globallity of a solution might be assumed it is hard to prove it. 4.3.2 Efficiency An important indicator on efficiency is the computing time a solver needs to solve an instance. This can be measured for example by the number of function evaluations, the number of iterations or simply in real time (given system specifications). In this research the number of iterations is used as a way to measure the efficiency of the solvers. Note that the number of iterations may not be a strong indicator, since for one solver a iteration may take longer than for the other. However, with this measure, we can make a distinction between easier or more difficult test cases. 23 Chapter 4: Solvers that the local solvers use different strategies since the found local optimum was not always the same for one test case. For test case 4, MINOS and SNOPT converged to an infeasible point, where CONOPT found the second best optimum. When one would want to find the global optimum of a bi-blending problem with several local optima, a local solver is not the appropriate optimization tool to use. Except for LGO, the other global optimizers were successful. In test case 4, BARON made some branching iterations; this was the only test case from this research which was ‘difficult’ for BARON in such a way that the global solution was not found (and verified) in the first iteration. OQNLP found all local optima and thus also found the global optimum to each test case. We conclude that BARON and OQNLP are the best solvers for solving the tested instances. 30 Chapter 5 Robustness 5.1 Introduction We would like to develop an algorithm that is able to identify solutions that have an -robustness with respect to the quadratic requirements. From practical considerations, one can define robustness R(x) of a design x∈X(2.4) with respect to Xas R(x) = max{R∈R+|(x+h)∈X, ∀h∈Rn,khk ≤ R}.(5.1) This -robustness concept is depicted in Figure 5.1. We want to find a solution to (bi-) blending while all solutions within a distance are feasible points as well. In the figure we can see that all points on the circle with a radius lie in the area g(x)≤0, and thus are feasible. The centre point of the circle is in this case an -robust point. g(x) = 0 X Figure 5.1: -robustness with respect to X. Solving the robustness of a point is a global optimization problem which can have local optima. We can measure the minimum distance of a point zto one quadratic constraint gpusing (5.2). 31 Chapter 5: Robustness rp(z) = min{pPn i=1(zi−xi)2} s.t. gp(x) = 0, p = 1, . . . , P (5.2) When we want to know the distance to the closest quadratic constraint in a global optimization problem we take the minimum over all distances rp,p= 1, . . . , P: R(x) = min p{rp(z)}, p = 1, . . . , P. (5.3) Notice that (5.2) is in general hard to solve; it may have several local optima and several KKT points that are not optimal [5]. 5.2 -Robustness with respect to (bi-)blending In this research all (bi-) blending problems were solved by standard solvers. It is difficult to combine the (bi-) blending problem while being robust simultaneously. In [4] a branch-and-bounds algorithm was designed to solve blending problems with a given -robustness. In our thesis however, the focus is on standard solvers which were really successful on non-robust (bi-) blending problems. That is why we try to reformulate the blending problem to be robust. In that case we would not have to change the optimization algorithm, but just the instances of the test case. In the next section this approach is explained in more detail. 5.3 Some (heuristic) approaches The idea of reformulating the blending problem is quite straight forward. The idea is to lay circles with a radius r=on the feasible side of a quadratic constraint and draw a new quadratic constraint through the centres of the circles. When we could formulate the problem like this, we could use our standard solvers to find an -robust optimum. In Figure 5.2 this approach is shown. The dotted lines are the ‘-robust constraints’. 32 Chapter 5: Robustness Figure 5.2: -robustness may be achieved by reformulation of quadratic constraints. The suggested approach would be a perfect method for solving the robustness problem. The approach leads us to the next question. How can we compute an -robust constraint? 5.3.1 Reformulation of quadratic constraints In Section 2.3 there was already a short introduction on quadratic functions. In this section we will give more insight on quadratic functions since we want to reformulate them. With this approach we have to distinguish definite and indefinite quadratic functions (see Figure 5.3). A quadratic function is definite when matrix Aof the quadratic form (5.4) has only positive (positive definite) or negative eigenvalues (negative definite). gp(x) = xTApx+bT px+dp(5.4) X X X X? X? X? X? Positive definite Negative definite Indefinite Figure 5.3: Definiteness of quadratic functions. 33 Chapter 5: Robustness We tried to reformulate the problem by adjusting d(see (5.4)) or by adjusting its eigenvalues. With the suggested method it appears impossible to lay equally sized circles on an indefinite constraint. This is because the limit of an indefinite function is always its eigenvector. In Figure 5.4 we sketch that laying equally sized circles on a constraint may be difficult for an indefinite function using the methods under investigation. 0 0 0 0 0 0 x1 x2 0 0 0 0 0 0.5 1 1.5 0 0.5 1 1.5 Figure 5.4: Circles with a different size on an indefinite function. That is why the focus in this chapter was on (negative) definite quadratic functions. At first, we tried to change a quadratic constraint by adjusting the constant term dfrom (5.4). We did this for a 2-dimensional instance g1with parameters: A1=1 0 0 3 , b1=0 0, d1=−1, d1=−0.5.(5.5) This gave us the result as shown in Figure 5.5 from which we can see that the found ellipse does not have equal distance from g1to g1for all points on the ellipse. Since distance ‘A’ is much bigger than distance ‘B’ this is not the appropriate reformulation. In Figure 5.6 the cross section of Figure 5.5 on g1is given in direction v2, with -robustness on v2. 34 Chapter 5: Robustness x1 x2 -1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 A B v2 v1 Figure 5.5: Reformulation by adjusting d. -1.2 0 g(x) KV2  Figure 5.6: Reformulation by adjusting d. [3] reformulates a quadratic function as (5.6). g(x)=(x−xc)TA(x−xc) + constant (5.6) in which xcis the stationary point of f(x) and constant is the function value of g(x) in the origin (g(xc)). One can derive that xcis given by (5.7) and constant is given by (5.8): xc=−1 2A−1b, (5.7) constant =d−1 4bTA−1b. (5.8) When we look at this notation we can see that only Adetermines the shape of the quadratic function, bdetermines (together with A) the origin and donly 35 Chapter 5: Robustness affects the function value (together with A) with the addition of a constant. The eigenvalues and eigenvectors determine the ‘shape’-characteristics of a quadratic function. The ‘eigenvalue decomposition’ of Matrix Ais written in (5.9): A=V DV T,(5.9) with Vrepresenting a matrix with the eigenvector vion Vi., and Da diagonal matrix with the eigenvalues λion the diagonal. When we decompose matrix A1from (5.5) we get: 1 0 0 3 =1 0 0 1  1 0 0 3  1 0 0 1 . From this decomposition we can conclude that on v1= (1, 0), λ1= 1. For v2 = (0, 1), λ2= 3. With this data we can adjust λ1and λ2, in such a way that the new λ’s can get us -robustness where the eigenvectors are intersecting g(x) = 0. In Figure 5.7 the cross section of Figure 5.5 on g1is given in direction v2, with -robustness on v2.kis the distance from xcto g(0) in the direction v2. -1.2 0 g(x) KV2  Figure 5.7: Cross section of g1on r2 By adjusting λwe can get the -robustness on r. This is derived in (5.10). k=q−g(xc) λ if pos.def: k=k−, if neg.def: k=k+ λ=−g(xc) k2  (5.10) When we do this for λ1and λ2from (5.5) we get the following robust matrix Dof adapted eigenvalues using = 0.05: D=1.1080 0 0 3.5959 . 36 Chapter 5: Robustness With (5.9) we obtain a matrix A1which is equal to Dsince Vis still the identity matrix. The ‘regular’ quadratic function g1is given by (5.5), while g1is given by: A1=1.1080 0 0 3.5959 , b1=0 0, d1=−1.(5.11) Figure 5.8 shows us that g1gives robust solutions on its eigenvectors. x1 x2 -1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 r 2 r 1 Figure 5.8: Robustness on the eigenvectors. The robust constraint seems to have an equal distance to the non-robust constraint as we desired. We can test this by solving the global optimization problem (5.2) with points zon g1. We can generate zby (5.12). We have to choose a vector vwhich is not an eigenvector of g(x). α=q−g(xc) VTAV z=αV (5.12) When we take VT= (1, 1) we obtain point zT= (0.4611, 0.4611) which lies on g1= 0. Substituting zin (5.2) together with g1gives a minimum robustness of 0.0490 when we solve this problem with Fmincon. This is smaller than the robustness we have obtained on the eigenvectors of 0.05. We can conclude that with this method we cannot generate an ellipse that is exactly -robust, but it comes close. 5.3.2 -robustness on X In the previous section we tried to reformulate quadratic constraints but we did not take the unit simplex into account. However, the feasible area of a blending instance is defined by X(2.4), which also includes the unit simplex. We can see what happens with the robustness on Xwhen we apply the eigenvalue approach from the previous section on test case 11 from Figure 2.6. In Figure 5.9 the 37 Chapter 5: Robustness robustness on the unit simplex in test case 11 is shown and it seems to be a good approximation. x1 x2 0 0.5 1 1.5 0 0.5 1 1.5 X X Figure 5.9: Robustness on test case 11. Another example shows us that we should include the unit simplex in the reformulation. This is shown in Figure 5.10. The difference between 1and 2is really big and not desirable. If we would want to solve this problem, maybe we should project the quadratic functions on the unit simplex beforehand using the ‘projection matrix’ from [2] . This is however not further researched in this thesis. Although the heuristic approach from this research does not guarantee solving the robustness problem with respect to (bi-) blending it can be used to generate starting points for a procedure. X1  X2  Figure 5.10: A big difference between 1and 2. 5.4 Conclusion Solving robustness on blending problems is a hard problem. Determining the robustness of a point is already a global optimization problem with several optima. 38 Chapter 5: Robustness When we add this problem to the (bi-) blending problem it gets really hard. In standard software it seems impossible to solve the problem directly. This is why we tried to reformulate the quadratic constraints so we could use the old problem formulation as used in Chapter 2 and 3. Changing the quadratic constraints by adjusting the constant term or the eigenvalue approach is not appropriate for reformulating indefinite quadratic functions. However, changing eigenvalues for (negative) definite functions seemed to work. However, we found a polynomial that was -robust on its eigenvectors, but not on the rest of its points. We neither took the unit simplex into account. In order to reformulate a blending problem you should consider the unit simplex at the same time since the feasible area is defined as X=Q∩S. This may be done using a projection matrix. 39 Appendix A: Testcases b4[3 ×1] =   0.0155 0.0515 −0.121   d4=−0.006 A5[3 ×3] =   20.605 −5.087 −10.9885 −5.087 32.003 −43.476 −10.9885 −43.476 −81.278   b5[3 ×1] =   0.1995 −0.097 126.7685   d5=−20.5063 A6[3 ×3] =   0.766 −0.1205 2.4735 −0.1205 0.528 1.9835 2.4735 1.9835 −7.822   b6[3 ×1] =   −2.432 −15.191 10.712   d6= 3.21125 A7[3 ×3] =   116.75 −3.09 168.553 −3.09 −67.424 515.114 168.553 515.114 −845.215   b7[3 ×1] =   −287.43 −645.926 354.537   46 Appendix A: Testcases d7= 115.0953 A8[3 ×3] =   1.0 3.0−0.5 3.0−5.0−3.5 −0.5−3.5−2.0   b8[3 ×1] =   0.832 0.832 0.832   d8= 0.968 A9[3 ×3] =   2.0−1.5 1.0 −1.5 1.0−1.0 1.0−1.0 3.0   b9[3 ×1] =   0.12 0.12 0.12   d9=−1.60 A10[3 ×3] =   4.0−1.5−1.5 −1.5 4.0−2.5 −1.5−2.5 4.0   b10[3 ×1] =   −0.026 −0.026 −0.026   d10 =−2.141 47 Appendix A: Testcases A11[3 ×3] =   −4.3 1.0−3.0 1.0 1.0−0.5 −3.0−0.5 2.5   b11[3 ×1] =   −0.193 −0.193 −0.193   d11 = 0.193 A12[3 ×3] =   4.0−1.0−2.0 −1.0 5.0−3.0 −2.0−3.0 4.0   b12[3 ×1] =   1.05 1.05 1.05   d12 =−2.1052 A13[3 ×3] =   −4.0 1.5 1.5 1.5−4.0 2.5 1.5 2.5−4.0   b13[3 ×1] =   −1.48 −1.48 −1.48   d13 = 2.36 A14[3 ×3] =   1.0 3.0−0.5 3.0−5.0−3.5 −0.5−3.5−2.0   48 Appendix A: Testcases b14[3 ×1] =   0.83 0.83 0.83   d14 =−0.59 A15[3 ×3] =   4.0−1.0 0.5 −1.0 1.0 3.5 0.5 3.5 0   b15[3 ×1] =   18.75 18.75 18.75   d15 =−19.5 A16[3 ×3] =   −2.0−0.5 0.5 −0.5 2.0−2.5 0.5−2.5 0   b16[3 ×1] =   10.89 10.89 10.89   d16 =−9.69 A17[3 ×3] =   −5.0 1.0−3 1 1 −0.5 −3−0.5 3   b17[3 ×1] =   6.45 6.45 6.45   49 Appendix A: Testcases d17 =−2.65 A18[3 ×3] =   5.0−1.0 1.0 −1.0 0 −2.5 1.0−2.5 1   b18[3 ×1] =   14.71 14.71 14.71   d18 =−14.88 A19[3 ×3] =   1.0 3.0−0.5 3.0−5.0−3.5 −0.5−3.5−2.0   b19[3 ×1] =   0.83 0.83 0.83   d19 =−0.95 A20[3 ×3] =   −4.0−3.0−0.5 −3.0−4.0−2.5 −0.5−2.5 3.0   b20[3 ×1] =   3.90 3.90 3.90   d20 =−0.69 50 Appendix A: Testcases A21[3 ×3] =   −2.0−0.5 0.5 −0.5 2.0−2.5 0.5−2.5 0   b21[3 ×1] =   2.08 2.08 2.08   d21 = 1.53 A22[3 ×3] =   −2.0−0.5 0.5 −0.5 2.0−2.5 0.5−2.5 0   b22[3 ×1] =   2.08 2.08 2.08   d22 =−0.87 A23[2 ×2] = −5.0 1.0 1.0−1.0 b23[2 ×1] = 3.0 0.5 d23 =−0.6 UniSpec1 Dimension = 7; Raw material cost = (114,115,107,127,115,106,108) Linear constraint: h1(x)=0.1493x1+ 0.6927x2+ 0.4643x3+ 0.7975x4+ 0.5967x5+ 0.6235x6+ 0.5284x7≥0.35 51 Appendix A: Testcases Quadratic constraints (gi(x) = xTAix+bT ix+di≤0; i= 1,2,3 ): A1[7×7] =           −1.473 8.215 −27.204 46.119 2.059 −11.929 −12.768 8.215 37.733346 5.127 95.691 34.954 20.165 19.445 −27.204 5.127 −21.743 36.843 −7.126 4.029 −4.152 46.119 95.691 36.843 189.643 93.359 52.904 54.802 2.059 34.954 −7.126 93.356 31.885 7.528 10.248 −11.929 20.165 4.029 52.904 7.528 11.951 10.964 −12.768 19.445 −4.152 54.802 10.248 10.964 7.197           b1[7 ×1] =           4.5675 34.7289 70.5707 −82.2761 29.3169 71.0818 63.7614           d1=−35 A2[7×7] =           1.35 −4.41 17.60 −92.45 2.74 −29.94 −14.05 −4.41 −39.13 −6.11 −126.38 −29.81 −63.42 −43.97 17.60 −6.11 15.45 −76.60 5.93 −44.05 −20.54 −92.45 −126.38 −76.60 −240.64 −117.46 −125.18 −114.98 2.74 −29.81 5.93 −117.46 −22.90 −47.37 −30.68 −29.94 −63.42 −44.05 −125.18 −47.37 −73.39 −73.99 −14.05 −43.97 −20.54 −114.98 −30.68 −73.99 −55.33           b2[7 ×1] =           −2.1232 −9.0403 −42.2072 190.5292 −9.9529 1.8162 5.1622           d2= 10 52 Appendix A: Testcases A3[7×7] =           −0.670 4.284 −12.837 23.708 1.677 −8.964 −4.859 4.284 21.380 −1.188 28.990 13.216 17.177 16.620 −12.837 −1.189 −21.376 9.841 −7.298 −10.043 −8.981 23.708 28.990 9.841 49.385 25.574 15.561 21.666 1.677 13.216 −7.298 25.574 8.419 4.149 6.595 −8.965 17.177 −10.043 15.561 4.149 1.090 6.292 −4.859 16.620 −8.981 21.666 6.594 6.292 5.906           b3[7 ×1] =           0.7097 −13.0982 27.5078 −49.1608 −7.3725 33.6731 11.3136           d3=−2 UniSpec5b Dimension = 7; Same raw material cost and similar quadratic requirements as UniSpec1 Quadratic constraints (gi(x) = xTAix+bT ix+di≤0; i= 4,5,6,7 ): A4=−A1;b4=−B1;d4= 45 A5=−A2;b5=−B2;d5=−21 A6[7 ×7] =           0.0−11.556 −1.114 14.690 −11.411 0.121 −0.150 −11.556 −3.316 −2.116 7.313 −8.800 19.897 9.051 −1.114 −2.116 4.728 16.250 −4.535 18.319 11.537 14.690 7.313 16.250 40.428 9.766 21.512 15.266 −11.412 −8.800 −4.535 9.766 −10.165 10.088 1.889 0.121 19.897 18.319 21.511 10.088 28.569 27.239 −0.150 9.051 11.537 15.266 1.889 27.239 19.965           53 Appendix A: Testcases b6[7 ×1] =           1.7278 23.5166 5.6724 −32.0798 19.0154 16.5074 7.31003           d6=−5 A7=A3;b7=B3;d7=−1 54 Appendix B GAMS code B.1 Rumcoke SETS I / i1, i2, i3 / P constraint / A, B /; alias (I,J); parameter L(I) lower bounds / i1 0, i2 0, i3 0 / U(I) upper bounds / i1 1, i2 1, i3 1 / D(P) constance D for function P / A -1, B 3.7/ C(I) cost per ingredient I / i1 0.1 i2 0.7 i3 4 /; table MATRIX(P,I,J) A Matrix for function P i1 i2 i3 A.i1 0 -16 0 A.i2 -16 0 0 A.i3 0 0 0 B.i1 10 0 2 B.i2 0 0 0 55