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