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 | |
---|