[37] | 1 | function massive |
---|
| 2 | ops{1} = sdpsettings('sos.cong',0,'sos.model',1,'verbose',0); |
---|
| 3 | ops{2} = sdpsettings('sos.cong',1,'sos.model',2,'verbose',0); |
---|
| 4 | ops{3} = sdpsettings('sos.cong',0,'sos.newton',0,'verbose',0,'sos.extlp',0); |
---|
| 5 | |
---|
| 6 | % Markus tentacles problem |
---|
| 7 | sdpvar x y a |
---|
| 8 | f = x^4 * y^2 + x^2 * y^4 - 3 * x^2 * y^2 + 1, k = 0; |
---|
| 9 | df = jacobian(f, [x y]); |
---|
| 10 | g = 1 - (df(1)^2 + df(2)^2) * (x^2 + y^2); |
---|
| 11 | if k >= 0 |
---|
| 12 | v = monolist([x; y], 2*k); |
---|
| 13 | coeffVec = sdpvar(length(v), 1); |
---|
| 14 | t = coeffVec' * v; |
---|
| 15 | constraints = set(sos(f - a - t * g)) + set(sos(t)); |
---|
| 16 | else |
---|
| 17 | coeffVec = []; |
---|
| 18 | constraints = set(sos(f - a)); |
---|
| 19 | end |
---|
| 20 | F = constraints; |
---|
| 21 | obj = -a; |
---|
| 22 | for i = 1:length(ops) |
---|
| 23 | i |
---|
| 24 | fail = regresstest(F,obj,ops{i},coeffVec); |
---|
| 25 | mbg_asserttolequal(fail,0); |
---|
| 26 | end |
---|
| 27 | |
---|
| 28 | |
---|
| 29 | |
---|
| 30 | |
---|
| 31 | function fail = regresstest(F,obj,ops,pv); |
---|
| 32 | |
---|
| 33 | if nargin==3 |
---|
| 34 | pv = []; |
---|
| 35 | end |
---|
| 36 | |
---|
| 37 | ops.sos.model = 1; |
---|
| 38 | solvesos(F,obj,ops,pv); |
---|
| 39 | obj1 = double(obj); |
---|
| 40 | p1s = checkset(F(find(is(F,'sos')))); |
---|
| 41 | p1e = checkset(F(find(~is(F,'sos')))); |
---|
| 42 | |
---|
| 43 | ops.sos.model = 2; |
---|
| 44 | solvesos(F,obj,ops,pv); |
---|
| 45 | obj2 = double(obj); |
---|
| 46 | p2s = checkset(F(find(is(F,'sos')))); |
---|
| 47 | p2e = checkset(F(find(~is(F,'sos')))); |
---|
| 48 | |
---|
| 49 | fail = 0; |
---|
| 50 | |
---|
| 51 | if abs(obj1-obj2) > 1e-4 |
---|
| 52 | fail = 1; |
---|
| 53 | end |
---|
| 54 | |
---|
| 55 | if any(p1s>1e-4) |
---|
| 56 | fail = 2; |
---|
| 57 | p1s |
---|
| 58 | end |
---|
| 59 | if any(p2s>1e-4) |
---|
| 60 | fail = 2; |
---|
| 61 | p2s |
---|
| 62 | end |
---|
| 63 | if any(p1e<-1e-4) |
---|
| 64 | fail = 2; |
---|
| 65 | p1e |
---|
| 66 | end |
---|
| 67 | if any(p2e<-1e-4) |
---|
| 68 | fail = 2; |
---|
| 69 | p2e |
---|
| 70 | end |
---|
| 71 | if fail==0 |
---|
| 72 | disp('Correct solution'); |
---|
| 73 | end |
---|