[37] | 1 | clc |
---|
| 2 | echo on |
---|
| 3 | %********************************************************* |
---|
| 4 | % |
---|
| 5 | % Moment relaxations |
---|
| 6 | % |
---|
| 7 | %********************************************************* |
---|
| 8 | % |
---|
| 9 | % This example shows how to solve some polynomial problems |
---|
| 10 | % using the theory of moment-relxations. |
---|
| 11 | pause |
---|
| 12 | clc |
---|
| 13 | |
---|
| 14 | % Define variables |
---|
| 15 | x1 = sdpvar(1,1); |
---|
| 16 | x2 = sdpvar(1,1); |
---|
| 17 | x3 = sdpvar(1,1); |
---|
| 18 | pause |
---|
| 19 | |
---|
| 20 | % Define a polynomial to be analyzed... |
---|
| 21 | p = -2*x1+x2-x3; |
---|
| 22 | pause |
---|
| 23 | |
---|
| 24 | % and a set of polynomial inequalities (concave constraint) |
---|
| 25 | F = set(x1*(4*x1-4*x2+4*x3-20)+x2*(2*x2-2*x3+9)+x3*(2*x3-13)+24>0); |
---|
| 26 | F = F + set(4-(x1+x2+x3)>0); |
---|
| 27 | F = F + set(6-(3*x2+x3)>0); |
---|
| 28 | F = F + set(x1>0); |
---|
| 29 | F = F + set(2-x1>0); |
---|
| 30 | F = F + set(x2>0); |
---|
| 31 | F = F + set(x3>0); |
---|
| 32 | F = F + set(3-x3>0); |
---|
| 33 | pause |
---|
| 34 | |
---|
| 35 | % A lower bound {min p(x) | F(x)>0} can be found using solvemoment |
---|
| 36 | solvemoment(F,p); |
---|
| 37 | pause |
---|
| 38 | clc |
---|
| 39 | % The lower bound can be be found by checking the |
---|
| 40 | % value of the RELAXED variables. When the moment relaxation |
---|
| 41 | % is calculated, x1 and x1^2 etc. are treated as independant |
---|
| 42 | % variables. To extract the values of the relaxed variables, |
---|
| 43 | % use the command relaxdouble |
---|
| 44 | relaxdouble(p) |
---|
| 45 | [double([x1;x2;x3;x1^2;x2^2;x3^2]) relaxdouble([x1;x2;x3;x1^2;x2^2;x3^2])] |
---|
| 46 | |
---|
| 47 | pause |
---|
| 48 | clc |
---|
| 49 | |
---|
| 50 | % Better lower bound can be obtained by using higher order |
---|
| 51 | % relaxations |
---|
| 52 | pause |
---|
| 53 | solvemoment(F,p,[],2); |
---|
| 54 | relaxdouble(p) |
---|
| 55 | |
---|
| 56 | pause |
---|
| 57 | clc |
---|
| 58 | |
---|
| 59 | % Another way to obtain better bounds is to add more valid constraints |
---|
| 60 | % One simple trick is to use the current lower bound by adding |
---|
| 61 | % the constraint p>double(p) |
---|
| 62 | pause |
---|
| 63 | solvemoment(F+set(p>double(p)),p,[],2); |
---|
| 64 | relaxdouble(p) |
---|
| 65 | |
---|
| 66 | pause |
---|
| 67 | clc |
---|
| 68 | % ...or square some linear constraints |
---|
| 69 | pause |
---|
| 70 | solvemoment(F+set(9>x3^2)+set(4>x1^2),p,[],2); |
---|
| 71 | relaxdouble(p) |
---|
| 72 | |
---|
| 73 | pause |
---|
| 74 | clc |
---|
| 75 | % or use an even higher relaxation |
---|
| 76 | pause |
---|
| 77 | solvemoment(F,p,[],4); |
---|
| 78 | relaxdouble(p) |
---|
| 79 | |
---|
| 80 | pause |
---|
| 81 | |
---|
| 82 | % The value -4 is known to be the global optimum, but our relaxed solution |
---|
| 83 | % is still not a valid solution. |
---|
| 84 | checkset(F) |
---|
| 85 | pause |
---|
| 86 | |
---|
| 87 | % YALMIP can try to recover globally optimal solutions, but to |
---|
| 88 | % do this, three output arguments are needed. The global solutions, |
---|
| 89 | % if sucessfully extracted, are returned in the second output. |
---|
| 90 | pause |
---|
| 91 | [solution,xoptimal] = solvemoment(F,p,[],4); |
---|
| 92 | |
---|
| 93 | % Use the first extracted global solution. |
---|
| 94 | setsdpvar([x1;x2;x3],xoptimal{1}); |
---|
| 95 | double(p) |
---|
| 96 | checkset(F) |
---|
| 97 | |
---|
| 98 | pause |
---|
| 99 | |
---|
| 100 | % In this example, there are only three variables and they |
---|
| 101 | % correpond to the variables in the output x. |
---|
| 102 | % |
---|
| 103 | % In principle, the code setsdpvar(recover([depends(p) depends(F)]),x{1}) |
---|
| 104 | % should work, but to be on the safe side in a more complex scenario, |
---|
| 105 | % a third output argument should be used to keep track of the variables. |
---|
| 106 | pause |
---|
| 107 | [solution,xoptimal,momentdata] = solvemoment(F,p,[],4); |
---|
| 108 | |
---|
| 109 | % The variables used in the relaxation are returned in |
---|
| 110 | % the fourth output, in the field 'x' |
---|
| 111 | momentdata |
---|
| 112 | pause |
---|
| 113 | |
---|
| 114 | setsdpvar(momentdata.x,xoptimal{1}); |
---|
| 115 | pause |
---|
| 116 | |
---|
| 117 | double(p) |
---|
| 118 | checkset(F) |
---|
| 119 | pause |
---|
| 120 | |
---|
| 121 | clc |
---|
| 122 | |
---|
| 123 | % Polynomial semidefinite constraints can be addressed also |
---|
| 124 | sdpvar x1 x2 |
---|
| 125 | p = -x1^2-x2^2; |
---|
| 126 | F = set([1-4*x1*x2 x1;x1 4-x1^2-x2^2]); |
---|
| 127 | pause |
---|
| 128 | [sol,xoptimal] = solvemoment(F,p,[],2); |
---|
| 129 | setsdpvar([x1;x2],xoptimal{1}); |
---|
| 130 | double(p) |
---|
| 131 | checkset(F) |
---|
| 132 | |
---|
| 133 | pause |
---|
| 134 | clc |
---|
| 135 | |
---|
| 136 | % Note that the moment relaxation solver can be called in |
---|
| 137 | % a fashion that is more consistent with the rest of the YALMIP framework. |
---|
| 138 | % |
---|
| 139 | % To do this, just use the solver tag 'moment'. |
---|
| 140 | pause |
---|
| 141 | solution = solvesdp(F,p,sdpsettings('solver','moment','moment.order',2)); |
---|
| 142 | setsdpvar(solution.momentdata.x,solution.xoptimal{1}); |
---|
| 143 | double(p) |
---|
| 144 | checkset(F) |
---|
| 145 | pause |
---|
| 146 | |
---|
| 147 | clc |
---|
| 148 | |
---|
| 149 | % The extraction of the global solution is numerically sensitive and |
---|
| 150 | % can easily lead to low accuracy, or even completly wrong solutions. |
---|
| 151 | % |
---|
| 152 | % Some tweaking can be done to improve upon this. |
---|
| 153 | % |
---|
| 154 | % Consider the following problem with known optima (0,2) and (0,-2) |
---|
| 155 | clear all |
---|
| 156 | yalmip('clear') |
---|
| 157 | x1 = sdpvar(1,1); |
---|
| 158 | x2 = sdpvar(1,1); |
---|
| 159 | obj = -x1^2-x2^2; |
---|
| 160 | g1 = 5-4*x1*x2-x1^2-x2^2; |
---|
| 161 | g2 = 4-16*x1*x2-2*x1^2-x2^2+4*x1^3*x2+4*x1*x2^3; |
---|
| 162 | F = set([g1 g2]); |
---|
| 163 | pause |
---|
| 164 | |
---|
| 165 | % We solve and extract solutions (requires a rather high order to find a solution) |
---|
| 166 | % Fo reasons that will become clear later, we also specify the tolerance |
---|
| 167 | % for the Gaussian elimination, used during the extraction process. |
---|
| 168 | pause |
---|
| 169 | [sol,xe] = solvemoment(F,obj,sdpsettings('moment.rceftol',1e-6),7); |
---|
| 170 | |
---|
| 171 | % The accuracy can easily become rather poor for this problem |
---|
| 172 | % NOTE : There is a randomization step in the extraction algorithm, |
---|
| 173 | % so these numbers may differ. |
---|
| 174 | [xe{1} xe{2}] |
---|
| 175 | pause |
---|
| 176 | |
---|
| 177 | % To improve the solutions, YALMIP can perform a number of Newton |
---|
| 178 | % iterations on a set of polynomial equalities that are solved using |
---|
| 179 | % numerical linear algebra techniques during the extraction algorithm. |
---|
| 180 | % (Note, these are not the original polynomials, but polynomials that |
---|
| 181 | % define the global solutions, given the moment matrices) |
---|
| 182 | pause |
---|
| 183 | [sol,xe] = solvemoment(F,obj,sdpsettings('moment.refine',5),7); |
---|
| 184 | [xe{1} xe{2}] |
---|
| 185 | pause |
---|
| 186 | |
---|
| 187 | % The main reason for the poor accuracy in this problem is actually |
---|
| 188 | % the bad conditioning during a Gaussian elimination. To improve |
---|
| 189 | % this situation, it is possible to let YALMIP select a tolerance using |
---|
| 190 | % some heuristics. |
---|
| 191 | % |
---|
| 192 | % This can be done by specifying the tolerance rceftol to -1. |
---|
| 193 | % In fact, this is the default setting, so we just run standard settings! |
---|
| 194 | pause |
---|
| 195 | [sol,xe] = solvemoment(F,obj,[],7); |
---|
| 196 | [xe{1} xe{2}] |
---|
| 197 | pause |
---|
| 198 | |
---|
| 199 | % For this problem we needed a rather high relaxation order to be able to |
---|
| 200 | % extract the solution. |
---|
| 201 | % |
---|
| 202 | % It is possible to tell YALMIP to try to extract a solution, even though |
---|
| 203 | % YALMIP belives this is impossible. |
---|
| 204 | % |
---|
| 205 | % To force YALMIP to try to extract global solutions, set the option |
---|
| 206 | % moment.extractrank to the desired number of solutions. |
---|
| 207 | % |
---|
| 208 | % Let us try to extract a solution from a 6th order relaxation. |
---|
| 209 | pause |
---|
| 210 | [sol,xe] = solvemoment(F,obj,sdpsettings('moment.extractrank',2),6); |
---|
| 211 | [xe{1} xe{2}] |
---|
| 212 | pause |
---|
| 213 | |
---|
| 214 | % Hmm, that's poor. No wonder, YALMIP did not expect it to be possible to |
---|
| 215 | % extract these solutions, so numerical problems are very likely during the |
---|
| 216 | % extraction process. Let us try to refine it. |
---|
| 217 | pause |
---|
| 218 | [sol,xe] = solvemoment(F,obj,sdpsettings('moment.extractrank',2,'moment.refine',5),6); |
---|
| 219 | [xe{1} xe{2}] |
---|
| 220 | pause |
---|
| 221 | |
---|
| 222 | % Pretty annoying to resolve the problem when you just want to change a |
---|
| 223 | % setting in the extraction algorithm, right? |
---|
| 224 | % |
---|
| 225 | % To avoid this, use 4 output arguments |
---|
| 226 | [sol,xe,momentdata] = solvemoment(F,obj,sdpsettings('moment.extractrank',2,'moment.refine',5),6); |
---|
| 227 | pause |
---|
| 228 | |
---|
| 229 | % The variable momentdata contains all information needed to extract |
---|
| 230 | % solutions. |
---|
| 231 | % |
---|
| 232 | % Hence, we can compare solutions for different settings easily. |
---|
| 233 | xe0 = extractsolution(momentdata,sdpsettings('moment.extractrank',2,'moment.refine',0)); |
---|
| 234 | xe1 = extractsolution(momentdata,sdpsettings('moment.extractrank',2,'moment.refine',2)); |
---|
| 235 | xe2 = extractsolution(momentdata,sdpsettings('moment.extractrank',2,'moment.refine',5)); |
---|
| 236 | [xe0{1} xe0{2} xe1{1} xe1{2} xe2{1} xe2{2}] |
---|
| 237 | pause |
---|
| 238 | |
---|
| 239 | echo off |
---|
| 240 | |
---|