[37] | 1 | function output = kktqp(interfacedata) |
---|
| 2 | %KKTQP Solver for indefinite QP problems using binary optmization |
---|
| 3 | |
---|
| 4 | % Author Johan Löfberg |
---|
| 5 | % $Id: kktqp.m,v 1.1 2006/04/04 19:01:38 joloef Exp $ |
---|
| 6 | |
---|
| 7 | % ******************************** |
---|
| 8 | %% INITIALIZE DIAGNOSTICS IN YALMIP |
---|
| 9 | % ******************************** |
---|
| 10 | K = interfacedata.K; |
---|
| 11 | b = interfacedata.F_struc(1+K.f:end,1); |
---|
| 12 | A = -interfacedata.F_struc(1+K.f:end,2:end); |
---|
| 13 | c = interfacedata.c; |
---|
| 14 | |
---|
| 15 | n = length(c); |
---|
| 16 | |
---|
| 17 | if ~all(isinf(interfacedata.lb)) |
---|
| 18 | if ~isempty(lb) |
---|
| 19 | A = [A;-eye(n)]; |
---|
| 20 | b = [b;-interfacedata.lb] |
---|
| 21 | end |
---|
| 22 | end |
---|
| 23 | if ~all(isinf(interfacedata.ub)) |
---|
| 24 | if ~isempty(ub) |
---|
| 25 | A = [A;eye(n)]; |
---|
| 26 | b = [b;interfacedata.ub] |
---|
| 27 | end |
---|
| 28 | end |
---|
| 29 | |
---|
| 30 | % Lazy... |
---|
| 31 | P = polytope(A,b); |
---|
| 32 | [B,L,U] = bounding_box(P); |
---|
| 33 | [A,b] = double(P); |
---|
| 34 | |
---|
| 35 | % Formulation here assumes maximization... |
---|
| 36 | Q = -2*interfacedata.Q; |
---|
| 37 | c = -interfacedata.c; |
---|
| 38 | |
---|
| 39 | x = sdpvar(n,1); |
---|
| 40 | bounds(x,L,U); |
---|
| 41 | |
---|
| 42 | y = sdpvar(length(b),1); % Duals |
---|
| 43 | dy = binvar(length(b),1); % indicater dual ==0 |
---|
| 44 | ds = binvar(length(b),1); % indicater slack==0 |
---|
| 45 | s = b-A*x; % slack |
---|
| 46 | |
---|
| 47 | % Derive bounds on primal slack |
---|
| 48 | [M,m] = derivebounds(s); |
---|
| 49 | |
---|
| 50 | % Let us try to derive bounds on the dauls |
---|
| 51 | % variables. Hmm... |
---|
| 52 | % Ugly, but let's just maximize/minimize linear |
---|
| 53 | % relaxations |
---|
| 54 | F = set(A'*y == Q*x + c) + set(s>0) + set(y>0);%KKT |
---|
| 55 | F = F + set(s < ds.*M); % Big M, we know upper bound on s |
---|
| 56 | F = F + set(dy+ds <= 1); % Complementary slackness |
---|
| 57 | F = F + set(0 <= sum(dy) <= n); |
---|
| 58 | |
---|
| 59 | % Find dis-joint constraints (silly way...) |
---|
| 60 | for i = 1:length(b) |
---|
| 61 | j = findrows(abs(A),abs(A(i,:))); |
---|
| 62 | if length(j)>1 |
---|
| 63 | S{i} = setdiff(j,i); |
---|
| 64 | else |
---|
| 65 | S{i} = []; |
---|
| 66 | end |
---|
| 67 | end |
---|
| 68 | |
---|
| 69 | [a1,a2,a3,model] = export(F,-y(i),sdpsettings('relax',1)); |
---|
| 70 | solvertime = clock; |
---|
| 71 | for i = 1:length(b) |
---|
| 72 | if isempty(S{i}) |
---|
| 73 | model_ = model; |
---|
| 74 | model_.c = model_.c*0; |
---|
| 75 | model_.c(n+i) = -1; |
---|
| 76 | sol = feval(model.solver.call,model_); |
---|
| 77 | else |
---|
| 78 | k = length(S{i}); |
---|
| 79 | Aeq = sparse(1:k,n+S{i},1,k,n+3*length(b)); |
---|
| 80 | beq = zeros(k,1); |
---|
| 81 | model_ = model; |
---|
| 82 | model_.F_struc = [beq Aeq;model.F_struc]; |
---|
| 83 | model_.K.f = model.K.f+k; |
---|
| 84 | model_.c = model_.c*0; |
---|
| 85 | model_.c(n+i) = -1; |
---|
| 86 | sol = feval(model.solver.call,model_); |
---|
| 87 | end |
---|
| 88 | if sol.problem == 0 |
---|
| 89 | My(i,1) = sol.Primal(n+i); |
---|
| 90 | else |
---|
| 91 | My(i,1) = 1e4; |
---|
| 92 | end |
---|
| 93 | end |
---|
| 94 | |
---|
| 95 | F = F + set(y <= dy.*My); |
---|
| 96 | |
---|
| 97 | obj = -0.5*(c'*x+b'*y); % ==cost in optimal points |
---|
| 98 | sol = solvesdp(F,obj); |
---|
| 99 | |
---|
| 100 | |
---|
| 101 | % ********************************** |
---|
| 102 | %% CREATE SOLUTION |
---|
| 103 | % ********************************** |
---|
| 104 | output.problem = sol.problem; |
---|
| 105 | output.Primal = double(x); |
---|
| 106 | output.Dual = []; |
---|
| 107 | output.Slack = []; |
---|
| 108 | output.infostr = yalmiperror(output.problem,'KKTQP'); |
---|
| 109 | output.solverinput = 0; |
---|
| 110 | if interfacedata.options.savesolveroutput |
---|
| 111 | output.solveroutput.solved_nodes = solved_nodes; |
---|
| 112 | output.solveroutput.lower = lower; |
---|
| 113 | output.solveroutput.upper = upper; |
---|
| 114 | else |
---|
| 115 | output.solveroutput =[]; |
---|
| 116 | end |
---|
| 117 | output.solvertime = etime(clock,solvertime); |
---|