[37] | 1 | function [Matrices,infeasible] = derive_bounds(Matrices,options) |
---|
| 2 | |
---|
| 3 | A = [Matrices.G -Matrices.E]; |
---|
| 4 | b = [Matrices.W]; |
---|
| 5 | Aeq = [Matrices.Aeq Matrices.Beq]; |
---|
| 6 | beq = [Matrices.beq]; |
---|
| 7 | lb = Matrices.lb; |
---|
| 8 | ub = Matrices.ub; |
---|
| 9 | |
---|
| 10 | oldlb = lb; |
---|
| 11 | oldub = ub; |
---|
| 12 | |
---|
| 13 | fixed = find(lb == ub); |
---|
| 14 | notfixed = setdiff(1:length(ub),fixed); |
---|
| 15 | |
---|
| 16 | if ~isempty(fixed) |
---|
| 17 | AA = A; |
---|
| 18 | bb = b; |
---|
| 19 | nn = lb(fixed); |
---|
| 20 | b = b-A(:,fixed)*lb(fixed); |
---|
| 21 | A(:,fixed) = []; |
---|
| 22 | |
---|
| 23 | lb = lb(notfixed); |
---|
| 24 | ub = ub(notfixed); |
---|
| 25 | end |
---|
| 26 | |
---|
| 27 | if ~isempty(fixed) & ~isempty(Aeq) |
---|
| 28 | beq = beq-Aeq(:,fixed)*nn; |
---|
| 29 | Aeq(:,fixed) = []; |
---|
| 30 | end |
---|
| 31 | |
---|
| 32 | n = size(A,2); |
---|
| 33 | In = eye(n); |
---|
| 34 | |
---|
| 35 | % Detect slacks and avoid maximizing them (they are unbounded) |
---|
| 36 | trivially_unbounded_above = zeros(n,1); |
---|
| 37 | trivially_unbounded_above(find(sum(A,1) == -sum(abs(A),1) & isinf(ub)')) = 1; |
---|
| 38 | trivially_unbounded_below = zeros(n,1); |
---|
| 39 | trivially_unbounded_below(find(sum(A,1) == sum(abs(A),1) & isinf(lb)')) = 1; |
---|
| 40 | if ~isempty(Aeq) |
---|
| 41 | trivially_unbounded_below(find(any(Matrices.Aeq,1))) = 0; |
---|
| 42 | trivially_unbounded_above(find(any(Matrices.Aeq,1))) = 0; |
---|
| 43 | end |
---|
| 44 | |
---|
| 45 | reset_to_minus_inf = find(isinf(ub(find(trivially_unbounded_below)))); |
---|
| 46 | reset_to_inf = find(isinf(ub(find(trivially_unbounded_above)))); |
---|
| 47 | |
---|
| 48 | lb=max(-1e7,lb); |
---|
| 49 | ub=min(1e7,ub); |
---|
| 50 | u = ub; |
---|
| 51 | l = lb; |
---|
| 52 | |
---|
| 53 | for i=1:n, |
---|
| 54 | |
---|
| 55 | if ~trivially_unbounded_above(i) |
---|
| 56 | f=-In(:,i); |
---|
| 57 | |
---|
| 58 | finite_lb = find(~isinf(l)); |
---|
| 59 | finite_ub = find(~isinf(u)); |
---|
| 60 | Abounds = [A;In(finite_ub,:);-In(finite_lb,:)]; |
---|
| 61 | bbounds = [b;u(finite_ub,:);-l(finite_lb,:)]; |
---|
| 62 | |
---|
| 63 | Abounds = [A;eye(n);-eye(n)]; |
---|
| 64 | bbounds = [b;u;-l]; |
---|
| 65 | |
---|
| 66 | [x,fval,lambda,exitflag,how]=mpt_solveLP(f(:)',Abounds,bbounds,[],[],[],options.mpt.lpsolver); |
---|
| 67 | if isequal(how,'infeasible') |
---|
| 68 | Matrices.lb = repmat(inf,length(oldlb),1); |
---|
| 69 | Matrices.ub = repmat(-inf,length(oldlb),1); |
---|
| 70 | infeasible = 1; |
---|
| 71 | return |
---|
| 72 | end |
---|
| 73 | u(i,1)=min(u(i,1),x(i)); |
---|
| 74 | end |
---|
| 75 | |
---|
| 76 | if ~trivially_unbounded_below(i) |
---|
| 77 | f=In(:,i); |
---|
| 78 | finite_lb = find(~isinf(l)); |
---|
| 79 | finite_ub = find(~isinf(u)); |
---|
| 80 | Abounds = [A;In(finite_ub,:);-In(finite_lb,:)]; |
---|
| 81 | bbounds = [b;u(finite_ub,:);-l(finite_lb,:)]; |
---|
| 82 | |
---|
| 83 | Abounds = [A;eye(n);-eye(n)]; |
---|
| 84 | bbounds = [b;u;-l]; |
---|
| 85 | |
---|
| 86 | [x,fval,lambda,exitflag,how]=mpt_solveLP(f(:)',Abounds,bbounds,[],[],[],options.mpt.lpsolver); |
---|
| 87 | if isequal(how,'infeasible') |
---|
| 88 | Matrices.lb = repmat(inf,length(oldlb),1); |
---|
| 89 | Matrices.ub = repmat(-inf,length(oldlb),1); |
---|
| 90 | infeasible = 1; |
---|
| 91 | return |
---|
| 92 | end |
---|
| 93 | l(i,1)=max(lb(i),x(i)); |
---|
| 94 | end |
---|
| 95 | |
---|
| 96 | end |
---|
| 97 | |
---|
| 98 | ll = oldlb; |
---|
| 99 | if ~isempty(notfixed) |
---|
| 100 | ll(notfixed) = l; |
---|
| 101 | end |
---|
| 102 | l = ll; |
---|
| 103 | uu = oldub; |
---|
| 104 | if ~isempty(notfixed) |
---|
| 105 | uu(notfixed) = u; |
---|
| 106 | end |
---|
| 107 | u = uu; |
---|
| 108 | infeasible = 0; |
---|
| 109 | |
---|
| 110 | % Now find bounds from equalities |
---|
| 111 | candidates = find(sum(Aeq | Aeq,2)==1); |
---|
| 112 | if length(candidates)>0 |
---|
| 113 | for i = candidates(:)' |
---|
| 114 | j = find(Aeq(i,:)); |
---|
| 115 | new_bound = beq(i)/Aeq(i,j); |
---|
| 116 | if any(u(j) < new_bound) | any(l(j)> new_bound) |
---|
| 117 | infeasible = 1; |
---|
| 118 | return |
---|
| 119 | else |
---|
| 120 | u(j)=new_bound; |
---|
| 121 | l(j)=new_bound; |
---|
| 122 | end |
---|
| 123 | end |
---|
| 124 | end |
---|
| 125 | |
---|
| 126 | trivially_unbounded_below = find(trivially_unbounded_below); |
---|
| 127 | trivially_unbounded_above = find(trivially_unbounded_above); |
---|
| 128 | lb(trivially_unbounded_below(reset_to_minus_inf)) = -inf; |
---|
| 129 | ub(trivially_unbounded_above(reset_to_inf)) = inf; |
---|
| 130 | Matrices.lb = lb; |
---|
| 131 | Matrices.ub = ub; |
---|