1 | function pout = propagatequadratics(p,upper,lower) |
---|
2 | |
---|
3 | pout = p; |
---|
4 | if p.bilinears~=0 |
---|
5 | F_struc = p.F_struc; |
---|
6 | |
---|
7 | p.F_struc = [-p.F_struc(1:p.K.f,:);p.F_struc]; |
---|
8 | p.K.f=2*p.K.f; |
---|
9 | % InequalityConstraintState = [ones(p.K.f,1);p.InequalityConstraintState]; |
---|
10 | InequalityConstraintState = [p.EqualityConstraintState;p.EqualityConstraintState;p.InequalityConstraintState]; |
---|
11 | |
---|
12 | if 0%upper < inf |
---|
13 | p.F_struc = [p.F_struc(1:p.K.f,:);upper-p.f -p.c';p.F_struc(1+p.K.f:end,:)]; |
---|
14 | p.K.l = p.K.l + 1; |
---|
15 | InequalityConstraintState = [1;InequalityConstraintState]; |
---|
16 | end |
---|
17 | |
---|
18 | if 0%~isnan(lower) |
---|
19 | p.F_struc = [-(p.lower-abs(p.lower)*0.01)+p.f p.c';p.F_struc]; |
---|
20 | InequalityConstraintState = [1;InequalityConstraintState]; |
---|
21 | p.K.l = p.K.l + 1; |
---|
22 | end |
---|
23 | |
---|
24 | if p.K.l+p.K.f>0 |
---|
25 | quadratic_variables = find(p.bilinears(:,2) == p.bilinears(:,3)); |
---|
26 | if ~isempty(quadratic_variables) |
---|
27 | quadratic_variables = p.bilinears(quadratic_variables,1); |
---|
28 | for i = 1:length(quadratic_variables) |
---|
29 | k = quadratic_variables(i); |
---|
30 | if p.lb(k) >= 0 & (p.lb(k) < p.ub(k)-1e-4) |
---|
31 | x = p.bilinears(p.bilinears(:,1)==k,2);% x^2 |
---|
32 | candidates = find((InequalityConstraintState==1) & p.F_struc(1:p.K.f+p.K.l,1+k))'; |
---|
33 | for j = candidates |
---|
34 | a = p.F_struc(j,2:end); |
---|
35 | aij = a(k); |
---|
36 | if aij > 0 |
---|
37 | indNEG = find(a < 0); |
---|
38 | indPOS = find(a > 0); |
---|
39 | LB = p.lb; |
---|
40 | UB = p.ub; |
---|
41 | LB(k) = 0; |
---|
42 | UB(k) = 0; |
---|
43 | a(k) = 0; |
---|
44 | newLB = (-p.F_struc(j,1)-a([indPOS(:);indNEG(:)])*[UB(indPOS);LB(indNEG)])/aij; |
---|
45 | p.lb(k) = max(p.lb(k),newLB); |
---|
46 | if p.lb(x)>0 |
---|
47 | p.lb(x) = max(p.lb(x),sqrt(max(0,newLB))); |
---|
48 | elseif p.ub(x)<0 |
---|
49 | p.ub(x) = min(p.ub(x),-sqrt(max(0,newLB))); |
---|
50 | end |
---|
51 | elseif aij < 0 |
---|
52 | indNEG = find(a < 0); |
---|
53 | indPOS = find(a > 0); |
---|
54 | LB = p.lb; |
---|
55 | UB = p.ub; |
---|
56 | LB(k) = 0; |
---|
57 | UB(k) = 0; |
---|
58 | a(k) = 0; |
---|
59 | newUB = (p.F_struc(j,1)+a([indPOS(:);indNEG(:)])*[UB(indPOS);LB(indNEG)])/(-aij); |
---|
60 | p.ub(k) = min(p.ub(k),newUB); |
---|
61 | p.ub(x) = min(p.ub(x),sqrt(max(0,newUB))); |
---|
62 | end |
---|
63 | end |
---|
64 | elseif p.lb(k)<0 |
---|
65 | % [p.lb(k) p.ub(k)] |
---|
66 | end |
---|
67 | end |
---|
68 | end |
---|
69 | end |
---|
70 | |
---|
71 | if p.K.l+p.K.f>0 |
---|
72 | bilinear_variables = find(p.bilinears(:,2) ~= p.bilinears(:,3)); |
---|
73 | if ~isempty(bilinear_variables) |
---|
74 | bilinear_variables = p.bilinears(bilinear_variables,1); |
---|
75 | for i = 1:length(bilinear_variables) |
---|
76 | k = bilinear_variables(i); |
---|
77 | if p.lb(k) >= -5000000000 & (p.lb(k) < p.ub(k)-1e-4) |
---|
78 | x = p.bilinears(p.bilinears(:,1)==k,2);% x^2 |
---|
79 | y = p.bilinears(p.bilinears(:,1)==k,3);% x^2 |
---|
80 | candidates = find((InequalityConstraintState==1) & p.F_struc(1:p.K.f+p.K.l,1+k))'; |
---|
81 | for j = candidates |
---|
82 | a = p.F_struc(j,2:end); |
---|
83 | aij = a(k); |
---|
84 | if aij > 0 |
---|
85 | indNEG = find(a < 0); |
---|
86 | indPOS = find(a > 0); |
---|
87 | LB = p.lb; |
---|
88 | UB = p.ub; |
---|
89 | LB(k) = 0; |
---|
90 | UB(k) = 0; |
---|
91 | a(k) = 0; |
---|
92 | newLB = (-p.F_struc(j,1)-a([indPOS(:);indNEG(:)])*[UB(indPOS);LB(indNEG)])/aij; |
---|
93 | p.lb(k) = max(p.lb(k),newLB); |
---|
94 | |
---|
95 | if p.lb(y) > 0 & p.ub(y)~=0 |
---|
96 | p.lb(x) = max(p.lb(x),newLB/p.ub(y)); |
---|
97 | end |
---|
98 | if p.lb(x) > 0 & p.ub(x)~=0 |
---|
99 | p.lb(y) = max(p.lb(y),newLB/p.ub(x)); |
---|
100 | end |
---|
101 | |
---|
102 | elseif aij < 0 |
---|
103 | indNEG = find(a < 0); |
---|
104 | indPOS = find(a > 0); |
---|
105 | LB = p.lb; |
---|
106 | UB = p.ub; |
---|
107 | LB(k) = 0; |
---|
108 | UB(k) = 0; |
---|
109 | a(k) = 0; |
---|
110 | newUB = (p.F_struc(j,1)+a([indPOS(:);indNEG(:)])*[UB(indPOS);LB(indNEG)])/(-aij); |
---|
111 | p.ub(k) = min(p.ub(k),newUB); |
---|
112 | if p.lb(y) > 0 & p.lb(y)~=0 |
---|
113 | p.ub(x) = min(p.ub(x),newUB/p.lb(y)); |
---|
114 | end |
---|
115 | if p.lb(x) > 0 & p.lb(x)~=0 |
---|
116 | p.ub(y) = min(p.ub(y),newUB/p.lb(x)); |
---|
117 | end |
---|
118 | end |
---|
119 | end |
---|
120 | elseif p.lb(k)<0 |
---|
121 | % [p.lb(k) p.ub(k)] |
---|
122 | end |
---|
123 | end |
---|
124 | end |
---|
125 | end |
---|
126 | end |
---|
127 | |
---|
128 | pout.lb = p.lb; |
---|
129 | pout.ub = p.ub; |
---|