[37] | 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 | |
---|