Mercurial > hg > octave-lyh

comparison scripts/optimization/qp.m @ 8280:5ee11a81688e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
author Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
date Tue, 28 Oct 2008 12:31:00 -0400
parents df9519e9990c
children cadc73247d65
comparison
equal deleted inserted replaced
8279:b3734f1cb592 8280:5ee11a81688e
129 n_in = 0; 129 n_in = 0;
130 if (nargin > 5) 130 if (nargin > 5)
131 if (! isempty (lb)) 131 if (! isempty (lb))
132 if (length(lb) != n) 132 if (length(lb) != n)
133 error ("qp: lower bound has incorrect length"); 133 error ("qp: lower bound has incorrect length");
134 else 134 elseif (isempty (ub))
135 Ain = [Ain; eye(n)]; 135 Ain = [Ain; eye(n)];
136 bin = [bin; lb]; 136 bin = [bin; lb];
137 endif 137 endif
138 endif 138 endif
139 139
140 if (! isempty (ub)) 140 if (! isempty (ub))
141 if (length (ub) != n) 141 if (length (ub) != n)
142 error ("qp: upper bound has incorrect length"); 142 error ("qp: upper bound has incorrect length");
143 else 143 elseif (isempty (lb))
144 Ain = [Ain; -eye(n)]; 144 Ain = [Ain; -eye(n)];
145 bin = [bin; -ub]; 145 bin = [bin; -ub];
146 endif 146 endif
147 endif
148
149 if (! isempty (lb) && ! isempty (ub))
150 rtol = sqrt (eps);
151 for i = 1:n
152 if (abs(lb (i) - ub(i)) < rtol*(1 + max (abs (lb(i) + ub(i)))))
153 ## These are actually an equality constraint
154 tmprow = zeros(1,n);
155 tmprow(i) = 1;
156 A = [A;tmprow];
157 b = [b; 0.5*(lb(i) + ub(i))];
158 n_eq = n_eq + 1;
159 else
160 tmprow = zeros(1,n);
161 tmprow(i) = 1;
162 Ain = [Ain; tmprow; -tmprow];
163 bin = [bin; lb(i); -ub(i)];
164 n_in = n_in + 2;
165 endif
166 endfor
147 endif 167 endif
148 endif 168 endif
149 169
150 ## Inequality constraints 170 ## Inequality constraints
151 if (nargin > 7) 171 if (nargin > 7)
154 error ("qp: inequality constraint matrix has incorrect column dimension"); 174 error ("qp: inequality constraint matrix has incorrect column dimension");
155 else 175 else
156 if (! isempty (A_lb)) 176 if (! isempty (A_lb))
157 if (length (A_lb) != dimA_in) 177 if (length (A_lb) != dimA_in)
158 error ("qp: inequality constraint matrix and lower bound vector inconsistent"); 178 error ("qp: inequality constraint matrix and lower bound vector inconsistent");
159 else 179 elseif (isempty (A_ub))
160 Ain = [Ain; A_in]; 180 Ain = [Ain; A_in];
161 bin = [bin; A_lb]; 181 bin = [bin; A_lb];
162 endif 182 endif
163 endif 183 endif
164 if (! isempty (A_ub)) 184 if (! isempty (A_ub))
165 if (length (A_ub) != dimA_in) 185 if (length (A_ub) != dimA_in)
166 error ("qp: inequality constraint matrix and upper bound vector inconsistent"); 186 error ("qp: inequality constraint matrix and upper bound vector inconsistent");
167 else 187 elseif (isempty (A_lb))
168 Ain = [Ain; -A_in]; 188 Ain = [Ain; -A_in];
169 bin = [bin; -A_ub]; 189 bin = [bin; -A_ub];
170 endif 190 endif
191 endif
192
193 if (! isempty (A_lb) && ! isempty (A_ub))
194 rtol = sqrt (eps);
195 for i=1:dimA_in
196 if (abs (A_lb(i) - A_ub(i)) < rtol*(1 + max (abs (A_lb(i) + A_ub(i)))))
197 ## These are actually an equality constraint
198 tmprow = A_in(i,:);
199 A = [A;tmprow];
200 b = [b; 0.5*(A_lb(i) + A_ub(i))];
201 n_eq = n_eq + 1;
202 else
203 tmprow = A_in(i,:);
204 Ain = [Ain; tmprow; -tmprow];
205 bin = [bin; A_lb(i); -A_ub(i)];
206 n_in = n_in + 2;
207 endif
208 endfor
171 endif 209 endif
172 endif 210 endif
173 endif 211 endif
174 212
175 ## Now we should have the following QP: 213 ## Now we should have the following QP:
193 rtol = sqrt (eps ("single")); 231 rtol = sqrt (eps ("single"));
194 else 232 else
195 rtol = sqrt (eps); 233 rtol = sqrt (eps);
196 endif 234 endif
197 235
198 eq_infeasible = (n_eq > 0 && norm (A*x0-b) > rtol*(1+norm (b))); 236 eq_infeasible = (n_eq > 0 && norm (A*x0-b) > rtol*(1+abs (b)));
199 in_infeasible = (n_in > 0 && any (Ain*x0-bin < -rtol*(1+norm (bin)))); 237 in_infeasible = (n_in > 0 && any (Ain*x0-bin < -rtol*(1+abs (bin))));
200 238
201 info = 0; 239 info = 0;
202 if (eq_infeasible || in_infeasible) 240 if (eq_infeasible || in_infeasible)
203 ## The initial guess is not feasible. 241 ## The initial guess is not feasible.
204 ## First define xbar that is feasible with respect to the equality 242 ## First define xbar that is feasible with respect to the equality
214 252
215 ## Check if xbar is feasible with respect to the inequality 253 ## Check if xbar is feasible with respect to the inequality
216 ## constraints also. 254 ## constraints also.
217 if (n_in > 0) 255 if (n_in > 0)
218 res = Ain * xbar - bin; 256 res = Ain * xbar - bin;
219 if (any (res < -rtol * (1 + norm (bin)))) 257 if (any (res < -rtol * (1 + abs (bin))))
220 ## xbar is not feasible with respect to the inequality 258 ## xbar is not feasible with respect to the inequality
221 ## constraints. Compute a step in the null space of the 259 ## constraints. Compute a step in the null space of the
222 ## equality constraints, by solving a QP. If the slack is 260 ## equality constraints, by solving a QP. If the slack is
223 ## small, we have a feasible initial guess. Otherwise, the 261 ## small, we have a feasible initial guess. Otherwise, the
224 ## problem is infeasible. 262 ## problem is infeasible.