1 | function [Matrices,infeasible] = mpt_reduce(Matrices) |
---|
2 | % Projects the whole mp(Q)LP problem on Aeq*U + Beq*x = beq |
---|
3 | % differs from mpt_project_on_equality in the sense that it |
---|
4 | % separates the integer/binary variables that have to be in |
---|
5 | % the basis |
---|
6 | % Aeq_cont*U_cont + Aeq_int*U_int + Beq*x = beq |
---|
7 | infeasible = 0; |
---|
8 | if length(Matrices.beq) > 0 |
---|
9 | |
---|
10 | [ii,jj,kk]=unique([Matrices.Aeq Matrices.Beq Matrices.beq],'rows'); |
---|
11 | |
---|
12 | integer_variables = union([Matrices.binary_variables Matrices.integer_variables]); |
---|
13 | cont_variables = setdiff(1:size(Matrices.Aeq,2),integer_variables); |
---|
14 | |
---|
15 | |
---|
16 | Matrices.Aeq = Matrices.Aeq(jj,:); |
---|
17 | |
---|
18 | Matrices.Aeq_cont = Matrices.Aeq(:,cont_variables); |
---|
19 | Matrices.Aeq_int = Matrices.Aeq(:,integer_variables); |
---|
20 | |
---|
21 | Matrices.Beq = Matrices.Beq(jj,:); |
---|
22 | Matrices.beq = Matrices.beq(jj,:); |
---|
23 | |
---|
24 | [Qh,Rh,e] = qr(full(Matrices.Aeq_cont),0); |
---|
25 | r = max(find(sum(abs(Rh),2)>1e-10)); |
---|
26 | |
---|
27 | % The dependent |
---|
28 | v1 = e(1:r); |
---|
29 | % The basis |
---|
30 | v2 = e(r+1:end); |
---|
31 | |
---|
32 | % H1u1+H2u2 = Mv + g |
---|
33 | Aeq1 = Matrices.Aeq_cont(:,v1); |
---|
34 | Aeq2 = Matrices.Aeq_cont(:,v2); |
---|
35 | |
---|
36 | Aeqtilde = [-Aeq1\Aeq2;eye(size(Aeq2,2))]; |
---|
37 | Beqtilde = [-Aeq1\Matrices.Beq;zeros(size(Aeq2,2),size(Matrices.Beq,2))]; |
---|
38 | beqtilde = [Aeq1\Matrices.beq;zeros(size(Aeq2,2),1)]; |
---|
39 | |
---|
40 | |
---|
41 | s = 1:size(Matrices.Aeq,2); |
---|
42 | p = zeros(1,length(s)); |
---|
43 | for i = 1:length(s) |
---|
44 | pi = find(s(i)==e); |
---|
45 | if ~isempty(pi) |
---|
46 | p(i) = pi; |
---|
47 | end |
---|
48 | end |
---|
49 | % This is what we would do in ML7.1 |
---|
50 | % [dummy,p] = ismember(1:size(Matrices.Aeq,2),e); |
---|
51 | |
---|
52 | S1 = Aeqtilde(p,:); |
---|
53 | S2 = Beqtilde(p,:); |
---|
54 | S3 = beqtilde(p,:); |
---|
55 | |
---|
56 | % New parameterization U = S1*z + S2*x + S3 |
---|
57 | M = Matrices; |
---|
58 | Matrices.G = M.G*S1; |
---|
59 | Matrices.E = M.E-M.G*S2; |
---|
60 | Matrices.W = M.W-M.G*S3; |
---|
61 | Matrices.nu = size(Matrices.G,2); |
---|
62 | |
---|
63 | if Matrices.qp |
---|
64 | Matrices.H = S1'*M.H*S1; |
---|
65 | Matrices.F = M.F*S1+S2'*M.H*S1; |
---|
66 | Matrices.Y = M.Y + S2'*M.H*S2+0.5*(M.F*S2+S2'*M.F'); |
---|
67 | Matrices.Cf = M.Cf*S1+S3'*M.H*S1; |
---|
68 | Matrices.Cc = M.Cc + M.Cf*S3; |
---|
69 | Matrices.Cx = M.Cx + S3'*M.F'+M.Cf*S2; |
---|
70 | else |
---|
71 | Matrices.H = M.H*S1; |
---|
72 | end |
---|
73 | |
---|
74 | removable = find(sum(abs([Matrices.G Matrices.E Matrices.G]),2)<1e-12); |
---|
75 | inconsistent = intersect(removable,find(Matrices.W<-1e-10)); |
---|
76 | if length(inconsistent)>0 |
---|
77 | infeasible = 1; |
---|
78 | return |
---|
79 | end |
---|
80 | |
---|
81 | if ~isempty(removable) |
---|
82 | Matrices.G(removable,:) = []; |
---|
83 | Matrices.E(removable,:) = []; |
---|
84 | Matrices.W(removable,:) = []; |
---|
85 | end |
---|
86 | |
---|
87 | % Keep the bounds for the new basis only |
---|
88 | Matrices.lb = [Matrices.lb(v2);Matrices.lb(end-size(Matrices.E,2)+1:end)]; |
---|
89 | Matrices.ub = [Matrices.ub(v2);Matrices.ub(end-size(Matrices.E,2)+1:end)]; |
---|
90 | |
---|
91 | % All equalities have been used |
---|
92 | Matrices.Aeq = []; |
---|
93 | Matrices.Beq = []; |
---|
94 | Matrices.beq = []; |
---|
95 | |
---|
96 | % This data is needed to recover original variables later |
---|
97 | if isempty(Matrices.getback) |
---|
98 | Matrices.getback.S1 = S1; |
---|
99 | Matrices.getback.S2 = S2; |
---|
100 | Matrices.getback.S3 = S3; |
---|
101 | else |
---|
102 | % This model has been reduced before, merge reductions |
---|
103 | oldgetback = Matrices.getback; |
---|
104 | Matrices.getback.S1 = oldgetback.S1*S1; |
---|
105 | Matrices.getback.S2 = oldgetback.S1*S2 + oldgetback.S2; |
---|
106 | Matrices.getback.S3 = oldgetback.S1*S3 + oldgetback.S3; |
---|
107 | end |
---|
108 | end |
---|
109 | |
---|