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