1 | function [sol, fval, exitflag, output] = ... |
---|
2 | optimize(funfcn, x0, lb, ub, A, b, Aeq, beq, nonlcon, strictness, options, varargin) |
---|
3 | %OPTIMIZE Optimize general constrained problems using FMINSEARCH |
---|
4 | % |
---|
5 | % Usage: |
---|
6 | % sol = OPTIMIZE(func, x0) |
---|
7 | % sol = OPTIMIZE(func, x0, lb, ub) |
---|
8 | % sol = OPTIMIZE(func, x0, lb, ub, A, b) |
---|
9 | % sol = OPTIMIZE(func, x0, lb, ub, A, b, Aeq, beq) |
---|
10 | % sol = OPTIMIZE(func, x0, lb, ub, A, b, Aeq, beq, nonlcon) |
---|
11 | % sol = OPTIMIZE(func, x0, lb, ub, A, b, Aeq, beq, nonlcon, strictness) |
---|
12 | % sol = OPTIMIZE(func, x0, lb, ub, A, b, Aeq, beq, nonlcon, strictness, options) |
---|
13 | % |
---|
14 | % [sol, fval] = OPTIMIZE(func, ...) |
---|
15 | % [sol, fval, exitflag] = OPTIMIZE(func, ...) |
---|
16 | % [sol, fval, exitflag, output] = OPTIMIZE(func, ...) |
---|
17 | % |
---|
18 | % INPUT ARGUMENTS: |
---|
19 | % |
---|
20 | % fun, x0, or options - see the help for FMINSEARCH |
---|
21 | % |
---|
22 | % lb - (OPTIONAL) lower bound vector or array, must have the same |
---|
23 | % size as x0. |
---|
24 | % |
---|
25 | % If no lower bounds exist for one of the variables, then |
---|
26 | % supply -inf for that variable. |
---|
27 | % |
---|
28 | % If no lower bounds exist at all, then [lb] may be left empty. |
---|
29 | % |
---|
30 | % Variables may be fixed in value by setting the corresponding |
---|
31 | % lower and upper bounds to exactly the same value. |
---|
32 | % |
---|
33 | % ub - (OPTIONAL) upper bound vector or array, must have the same |
---|
34 | % size as x0. |
---|
35 | % |
---|
36 | % If no upper bounds exist for one of the variables, then |
---|
37 | % supply +inf for that variable. |
---|
38 | % |
---|
39 | % If no upper bounds at all, then [ub] may be left empty. |
---|
40 | % |
---|
41 | % Variables may be fixed in value by setting the corresponding |
---|
42 | % lower and upper bounds to exactly the same value. |
---|
43 | % |
---|
44 | % A, b - (OPTIONAL) Linear inequality constraint array and right |
---|
45 | % hand side vector. (Note: these constraints were chosen to |
---|
46 | % be consistent with those of fmincon.) |
---|
47 | % |
---|
48 | % This linear constraint forces the solution vector [x] to |
---|
49 | % satisfy |
---|
50 | % A*x <= b |
---|
51 | % |
---|
52 | % Note that in case [x] is a matrix (this is true when [x0] is |
---|
53 | % a matrix), the argument [b] must have corresponding size |
---|
54 | % [size(A,1) x size(x0,2)], since the same equation is used to |
---|
55 | % evaluate this constraint. |
---|
56 | % |
---|
57 | % Aeq, beq - (OPTIONAL) Linear equality constraint array and right |
---|
58 | % hand side vector. (Note: these constraints were chosen to |
---|
59 | % be consistent with those of fmincon.) |
---|
60 | % |
---|
61 | % This linear constraint forces the solution vector [x] to |
---|
62 | % satisfy |
---|
63 | % |
---|
64 | % Aeq*x == beq |
---|
65 | % |
---|
66 | % Note that in case [x] is a matrix (this is true when [x0] is |
---|
67 | % a matrix), the argument [beq] must have corresponding size |
---|
68 | % [size(Aeq,1) x size(x0,2)], since the same equation is used to |
---|
69 | % evaluate this constraint. |
---|
70 | % |
---|
71 | % nonlcon - (OPTIONAL) function handle to general nonlinear constraints, |
---|
72 | % inequality and/or equality constraints. |
---|
73 | % |
---|
74 | % [nonlcon] must return two vectors, [c] and [ceq], containing the |
---|
75 | % values for the nonlinear inequality constraints [c] and |
---|
76 | % those for the nonlinear equality constraints [ceq] at [x]. (Note: |
---|
77 | % these constraints were chosen to be consistent with those of |
---|
78 | % fmincon.) |
---|
79 | % |
---|
80 | % These constraints force the solution to satisfy |
---|
81 | % |
---|
82 | % ceq(x) = 0 |
---|
83 | % c(x) <= 0, |
---|
84 | % |
---|
85 | % where [c(x)] and [ceq(x)] are general non-linear functions of [x]. |
---|
86 | % |
---|
87 | % strictness - (OPTIONAL) By default, OPTIMIZE will assume the objective |
---|
88 | % (and constraint) function(s) can be evaluated at ANY point in |
---|
89 | % RN-space; the initial estimate does not have to lie in the |
---|
90 | % feasible region, and intermediate solutions are also allowed to step |
---|
91 | % outside this area. If your function does not permit such behavior, |
---|
92 | % set this argument to 'strict'. With 'strict' enabled, the linear |
---|
93 | % constraints will be satisfied strictly, while the nonlinear |
---|
94 | % constraints will be satisfied within options.TolCon. |
---|
95 | % |
---|
96 | % If this is also not permissible, use 'superstrict' - then all |
---|
97 | % nonlinear constraints are also satisfied AT ALL TIMES, and the |
---|
98 | % objective function is NEVER evaluated outside the feasible area. |
---|
99 | % |
---|
100 | % When using 'strict' or 'superstrict', the initial estimate [x0] |
---|
101 | % MUST be feasible. If it is not feasible, an error is produced |
---|
102 | % before the objective function is ever evaluated. |
---|
103 | % |
---|
104 | % |
---|
105 | % OUTPUT ARGUMENTS: |
---|
106 | % |
---|
107 | % sol, fval - the solution vector and the corresponding function value, |
---|
108 | % respectively. |
---|
109 | % |
---|
110 | % exitflag - (See also the help on FMINSEARCH) A flag that specifies the |
---|
111 | % reason the algorithm terminated. FMINSEARCH uses only the values |
---|
112 | % |
---|
113 | % 1 fminsearch converged to a solution x |
---|
114 | % 0 Max. # of function evaluations or iterations exceeded |
---|
115 | % -1 Algorithm was terminated by the output function. |
---|
116 | % |
---|
117 | % Since OPTIMIZE handles constrained problems, the following two |
---|
118 | % values were added: |
---|
119 | % |
---|
120 | % 2 All elements in [lb] and [ub] were equal - nothing done |
---|
121 | % -2 Problem is infeasible after the optimization (Some or |
---|
122 | % any of the constraints are violated at the final |
---|
123 | % solution). |
---|
124 | % |
---|
125 | % output - (See also the help on FMINSEARCH) A structure that contains |
---|
126 | % additional details on the optimization. FMINSEARCH returns |
---|
127 | % |
---|
128 | % output.algorithm Algorithm used |
---|
129 | % output.funcCount Number of function evaluations |
---|
130 | % output.iterations Number of iterations |
---|
131 | % output.message Exit message |
---|
132 | % |
---|
133 | % Since OPTIMIZE handles constrained problems, the following |
---|
134 | % fields were added: |
---|
135 | % |
---|
136 | % output.constrviolation.lin_ineq |
---|
137 | % output.constrviolation.lin_eq |
---|
138 | % output.constrviolation.nonlin_ineq |
---|
139 | % output.constrviolation.nonlin_ineq |
---|
140 | % |
---|
141 | % All these fields contain a [M x 2]-cell array. The fist column |
---|
142 | % contains a logical index to the constraints, which is true if the |
---|
143 | % constraint was violated, false if it was satisfied. The second |
---|
144 | % column contains the amount of constraint violation. This amount is |
---|
145 | % equal to zero if the constraint was satisfied within |
---|
146 | % options.TolCon. |
---|
147 | % |
---|
148 | % |
---|
149 | % Notes: |
---|
150 | % |
---|
151 | % If options is supplied, then TolX will apply to the transformed |
---|
152 | % variables. All other FMINSEARCH parameters should be unaffected. |
---|
153 | % |
---|
154 | % Variables which are constrained by both a lower and an upper |
---|
155 | % bound will use a sin() transformation. Those constrained by |
---|
156 | % only a lower or an upper bound will use a quadratic |
---|
157 | % transformation, and unconstrained variables will be left alone. |
---|
158 | % |
---|
159 | % Variables may be fixed by setting their respective bounds equal. |
---|
160 | % In this case, the problem will be reduced in size for FMINSEARCH. |
---|
161 | % |
---|
162 | % If your problem has an EXCLUSIVE (strict) bound constraints which |
---|
163 | % will not permit evaluation at the bound itself, then you must |
---|
164 | % provide a slightly offset bound. An example of this is a function |
---|
165 | % which contains the log of one of its parameters. If you constrain |
---|
166 | % the variable to have a lower bound of zero, then OPTIMIZE may |
---|
167 | % try to evaluate the function exactly at zero. |
---|
168 | % |
---|
169 | % EXAMPLES: |
---|
170 | % |
---|
171 | % rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2; |
---|
172 | % |
---|
173 | % <<Fully unconstrained problem>> |
---|
174 | % |
---|
175 | % optimize(rosen, [3 3]) |
---|
176 | % ans = |
---|
177 | % 1.0000 1.0000 |
---|
178 | % |
---|
179 | % |
---|
180 | % <<lower bound constrained>> |
---|
181 | % |
---|
182 | % optimize(rosen,[3 3],[2 2],[]) |
---|
183 | % ans = |
---|
184 | % 2.0000 4.0000 |
---|
185 | % |
---|
186 | % |
---|
187 | % <<x(2) fixed at 3>> |
---|
188 | % |
---|
189 | % optimize(rosen,[3 3],[-inf 3],[inf,3]) |
---|
190 | % ans = |
---|
191 | % 1.7314 3.0000 |
---|
192 | % |
---|
193 | % |
---|
194 | % <<simple linear inequality: x(1) + x(2) <= 1>> |
---|
195 | % |
---|
196 | % optimize(rosen,[0 0],[],[],[1 1], 1) |
---|
197 | % |
---|
198 | % ans = |
---|
199 | % 0.6187 0.3813 |
---|
200 | % |
---|
201 | % |
---|
202 | % <<nonlinear inequality: sqrt(x(1)^2 + x(2)^2) <= 1>> |
---|
203 | % <<nonlinear equality : x(1)^2 + x(2)^3 = 0.5>> |
---|
204 | % |
---|
205 | % execute this m-file: |
---|
206 | % |
---|
207 | % function test_optimize |
---|
208 | % rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2; |
---|
209 | % |
---|
210 | % options = optimset('TolFun', 1e-8, 'TolX', 1e-8); |
---|
211 | % |
---|
212 | % optimize(rosen, [3 3], [],[],[],[],[],[],... |
---|
213 | % @nonlcon, [], options) |
---|
214 | % |
---|
215 | % end |
---|
216 | % function [c, ceq] = nonlcon(x) |
---|
217 | % c = norm(x) - 1; |
---|
218 | % ceq = x(1)^2 + x(2)^3 - 0.5; |
---|
219 | % end |
---|
220 | % |
---|
221 | % ans = |
---|
222 | % 0.6513 0.4233 |
---|
223 | % |
---|
224 | % |
---|
225 | % |
---|
226 | % Of course, any combination of the above constraints is |
---|
227 | % also possible. |
---|
228 | % |
---|
229 | % |
---|
230 | % See also: fminsearch, fminsearchcon, fminsearchbnd, fmincon. |
---|
231 | |
---|
232 | |
---|
233 | % Author : Rody P.S. Oldenhuis |
---|
234 | % Delft University of Technology |
---|
235 | % E-mail : oldenhuis@dds.nl |
---|
236 | % |
---|
237 | % FMINSEARCHBND, FMINSEARCHCON and most of the |
---|
238 | % help for OPTIMIZE witten by |
---|
239 | % : John D'Errico |
---|
240 | % E-mail : woodchips@rochester.rr.com |
---|
241 | % |
---|
242 | % Last edited: 29/May/2009 |
---|
243 | % Uploaded : 29/May/2009 |
---|
244 | % |
---|
245 | % [ Please report bugs to oldenhuis@dds.nl ] |
---|
246 | |
---|
247 | %% initialize |
---|
248 | |
---|
249 | % process input |
---|
250 | error(nargchk(2, inf, nargin)) |
---|
251 | % FUTURE WORK: no x0 given (nargin < 2) means optimize globally. |
---|
252 | % for now, just produce an error in that case |
---|
253 | if (nargin < 11), options = optimset; end |
---|
254 | if (nargin < 10), strictness = 'loose'; end |
---|
255 | if (nargin < 9), nonlcon = []; end |
---|
256 | if (nargin < 8), beq = []; end |
---|
257 | if (nargin < 7), Aeq = []; end |
---|
258 | if (nargin < 6), b = []; end |
---|
259 | if (nargin < 5), A = []; end |
---|
260 | if (nargin < 4), lb = []; end |
---|
261 | if (nargin < 3), ub = []; end |
---|
262 | |
---|
263 | % get tolerance on constraints |
---|
264 | tolCon = optimget(options, 'TolCon', 1e-6); |
---|
265 | |
---|
266 | % check for an output function. If there is any, |
---|
267 | % use wrapper function to call it with untransformed variable |
---|
268 | if ~isempty(options.OutputFcn) |
---|
269 | OutputFcn = options.OutputFcn; |
---|
270 | options.OutputFcn = @OutputFcn_wrapper; |
---|
271 | end |
---|
272 | |
---|
273 | % define number of dimensions |
---|
274 | N = numel(x0); |
---|
275 | |
---|
276 | % adjust bounds when they are empty |
---|
277 | if isempty(lb), lb = -inf(size(x0)); end |
---|
278 | if isempty(ub), ub = +inf(size(x0)); end |
---|
279 | |
---|
280 | % check the user-provided input with nested function check_input |
---|
281 | check_input; |
---|
282 | |
---|
283 | % resize and reshape all input to be of completely predictable size |
---|
284 | new_x = x0; ub = ub(:); |
---|
285 | x0 = x0(:); lb = lb(:); |
---|
286 | |
---|
287 | % replicate lb or ub when they are scalars, and x0 is not |
---|
288 | if isscalar(lb) && (N ~= 1), lb = repmat(lb, size(x0)); end |
---|
289 | if isscalar(ub) && (N ~= 1), ub = repmat(ub, size(x0)); end |
---|
290 | |
---|
291 | % determine the type of bounds |
---|
292 | nf_lb = ~isfinite(lb); nf_ub = ~isfinite(ub); |
---|
293 | lb_only = ~nf_lb & nf_ub; ub_only = nf_lb & ~nf_ub; |
---|
294 | unconst = nf_lb & nf_ub; fix_var = lb == ub; |
---|
295 | lb_ub = ~nf_lb & ~nf_ub & ~fix_var; |
---|
296 | |
---|
297 | % if all variables are fixed, simply return |
---|
298 | if length(lb(fix_var)) == N |
---|
299 | sol = reshape(lb,size(newx)); output.iterations = 0; |
---|
300 | output.funcCount = 0; fval = funfcn(sol); |
---|
301 | output.message = 'Lower and upper bound were set equal - nothing to do. '; |
---|
302 | exitflag = 2; [output, exitflag] = finalize(lb, output, exitflag); |
---|
303 | if exitflag ~= -2 |
---|
304 | output.message = sprintf(... |
---|
305 | '%s\nFortunately, the solution is feasible using OPTIONS.TolCon of %1.6f.',... |
---|
306 | output.message, tolCon); |
---|
307 | end, return; |
---|
308 | end |
---|
309 | |
---|
310 | % force the initial estimate inside the given bounds |
---|
311 | x0(x0 < lb) = lb(x0 < lb); x0(x0 > ub) = ub(x0 > ub); |
---|
312 | |
---|
313 | % transform initial estimate to its unconstrained counterpart |
---|
314 | xin = x0; % fixed and unconstrained variables |
---|
315 | xin(lb_only) = sqrt(x0(lb_only) - lb(lb_only)); % lower bounds only |
---|
316 | xin(ub_only) = sqrt(ub(ub_only) - x0(ub_only)); % upper bounds only |
---|
317 | xin(lb_ub) = real(asin( 2*(x0(lb_ub) - lb(lb_ub))./ ... |
---|
318 | (ub(lb_ub) - lb(lb_ub)) - 1)); % both upper and lower bounds |
---|
319 | |
---|
320 | % remove fixed variables from the initial estimate |
---|
321 | xin(fix_var) = []; |
---|
322 | |
---|
323 | % some constant matrices to speed things up (slightly) |
---|
324 | Nzero = zeros(N, 1); Np1zero = zeros(N, N+1); exp200 = exp(200); |
---|
325 | |
---|
326 | % define the transformed & penalized function |
---|
327 | funfncT = @(x) funfncP(X(x), varargin{:}); |
---|
328 | |
---|
329 | % optimize the problem with FMINSEARCH |
---|
330 | [presol, fval, exitflag, output] = fminsearch(funfncT, xin, options, varargin{:}); |
---|
331 | |
---|
332 | % and transform everything back to original (bounded) variables |
---|
333 | sol = new_x; sol(:) = X(presol); % with the same size as the original x0 |
---|
334 | |
---|
335 | % append constraint violations to the output structure, and change the |
---|
336 | % exitflag accordingly |
---|
337 | [output, exitflag] = finalize(sol, output, exitflag); |
---|
338 | |
---|
339 | %% nested functions (the actual work) |
---|
340 | |
---|
341 | % check user provided input |
---|
342 | function check_input |
---|
343 | |
---|
344 | % dimensions & weird input |
---|
345 | if (numel(lb) ~= N && ~isscalar(lb)) || (numel(ub) ~= N && ~isscalar(ub)) |
---|
346 | error('optimize:lb_ub_incompatible_size',... |
---|
347 | 'Size of either [lb] or [ub] incompatible with size of [x0].') |
---|
348 | end |
---|
349 | if ~isempty(A) && isempty(b) |
---|
350 | warning('optimize:Aeq_but_not_beq', ... |
---|
351 | ['I received the matrix [A], but you omitted the corresponding vector [b].',... |
---|
352 | '\nI`ll assume a zero-vector for [b]...']) |
---|
353 | b = zeros(size(A,1), size(x0,2)); |
---|
354 | end |
---|
355 | if ~isempty(Aeq) && isempty(beq) |
---|
356 | warning('optimize:Aeq_but_not_beq', ... |
---|
357 | ['I received the matrix [Aeq], but you omitted the corresponding vector [beq].',... |
---|
358 | '\nI`ll assume a zero-vector for [beq]...']) |
---|
359 | beq = zeros(size(Aeq,1), size(x0,2)); |
---|
360 | end |
---|
361 | if isempty(Aeq) && ~isempty(beq) |
---|
362 | warning('optimize:beq_but_not_Aeq', ... |
---|
363 | ['I received the vector [beq], but you omitted the corresponding matrix [Aeq].',... |
---|
364 | '\nI`ll ignore the given [beq]...']) |
---|
365 | beq = []; |
---|
366 | end |
---|
367 | if isempty(A) && ~isempty(b) |
---|
368 | warning('optimize:b_but_not_A', ... |
---|
369 | ['I received the vector [b], but you omitted the corresponding matrix [A].',... |
---|
370 | '\nI`ll ignore the given [b]...']) |
---|
371 | b = []; |
---|
372 | end |
---|
373 | if ~isempty(A) && ~isempty(b) && size(b,1)~=size(A,1) |
---|
374 | error('optimize:b_incompatible_with_A',... |
---|
375 | 'The size of [b] is incompatible with that of [A].') |
---|
376 | end |
---|
377 | if ~isempty(Aeq) && ~isempty(beq) && size(beq,1)~=size(Aeq,1) |
---|
378 | error('optimize:b_incompatible_with_A',... |
---|
379 | 'The size of [beq] is incompatible with that of [Aeq].') |
---|
380 | end |
---|
381 | if ~isvector(x0) && ~isempty(A) && (size(A,2) ~= size(x0,1)) |
---|
382 | error('optimize:A_incompatible_size',... |
---|
383 | 'Linear constraint matrix [A] has incompatible size for given [x0].') |
---|
384 | end |
---|
385 | if ~isvector(x0) && ~isempty(Aeq) && (size(Aeq,2) ~= size(x0,1)) |
---|
386 | error('optimize:Aeq_incompatible_size',... |
---|
387 | 'Linear constraint matrix [Aeq] has incompatible size for given [x0].') |
---|
388 | end |
---|
389 | if ~isempty(b) && size(b,2)~=size(x0,2) |
---|
390 | error('optimize:x0_vector_but_not_b',... |
---|
391 | 'Given linear constraint vector [b] has incompatible size with given [x0].') |
---|
392 | end |
---|
393 | if ~isempty(beq) && size(beq,2)~=size(x0,2) |
---|
394 | error('optimize:x0_vector_but_not_beq',... |
---|
395 | 'Given linear constraint vector [beq] has incompatible size with given [x0].') |
---|
396 | end |
---|
397 | |
---|
398 | % functions are not function handles |
---|
399 | if ~isa(funfcn, 'function_handle') |
---|
400 | error('optimize:func_not_a_function',... |
---|
401 | 'Objective function must be given as a function handle.') |
---|
402 | end |
---|
403 | if ~isempty(nonlcon) && ~isa(nonlcon, 'function_handle') |
---|
404 | error('optimize:nonlcon_not_a_function', ... |
---|
405 | 'non-linear constraint function must be a function handle.') |
---|
406 | end |
---|
407 | |
---|
408 | % test the feasibility of the initial solution (when strict or |
---|
409 | % superstrict behavior has been enabled) |
---|
410 | if strcmpi(strictness, 'strict') || strcmpi(strictness, 'superstrict') |
---|
411 | superstrict = strcmpi(strictness, 'superstrict'); |
---|
412 | if ~isempty(A) && any(any(A*x0 > b)) |
---|
413 | error('optimize:x0_doesnt_satisfy_linear_ineq', ... |
---|
414 | ['Initial estimate does not satisfy linear inequality.', ... |
---|
415 | '\nPlease provide an initial estimate inside the feasible region.']); |
---|
416 | end |
---|
417 | if ~isempty(Aeq) && any(any(Aeq*x0 ~= beq)) |
---|
418 | error('optimize:x0_doesnt_satisfy_linear_eq', ... |
---|
419 | ['Initial estimate does not satisfy linear equality.', ... |
---|
420 | '\nPlease provide an initial estimate inside the feasible region.']); |
---|
421 | end |
---|
422 | if ~isempty(nonlcon) |
---|
423 | [c, ceq] = nonlcon(x0); %#ok |
---|
424 | if ~isempty(c) && any(c(:) > ~superstrict*tolCon) |
---|
425 | error('optimize:x0_doesnt_satisfy_nonlinear_ineq', ... |
---|
426 | ['Initial estimate does not satisfy nonlinear inequality.', ... |
---|
427 | '\nPlease provide an initial estimate inside the feasible region.']); |
---|
428 | end |
---|
429 | if ~isempty(ceq) && any(abs(ceq(:)) >= ~superstrict*tolCon) |
---|
430 | error('optimize:x0_doesnt_satisfy_nonlinear_eq', ... |
---|
431 | ['Initial estimate does not satisfy nonlinear equality.', ... |
---|
432 | '\nPlease provide an initial estimate inside the feasible region.']); |
---|
433 | end |
---|
434 | end |
---|
435 | end |
---|
436 | |
---|
437 | end % check_input |
---|
438 | |
---|
439 | % create transformed variable to conform to upper and lower bounds |
---|
440 | function xx = X(x) |
---|
441 | if (size(x, 2) == 1), xx = Nzero; else xx = Np1zero; end |
---|
442 | for i = 1:size(x, 2) |
---|
443 | xx(lb_only, i) = lb(lb_only) + x(lb_only(~fix_var), i).^2; |
---|
444 | xx(ub_only, i) = ub(ub_only) - x(ub_only(~fix_var), i).^2; |
---|
445 | xx(fix_var, i) = lb(fix_var); |
---|
446 | xx(unconst, i) = x(unconst(~fix_var), i); |
---|
447 | xx(lb_ub, i) = lb(lb_ub) + ... |
---|
448 | (ub(lb_ub) - lb(lb_ub)) .* (sin(x(lb_ub(~fix_var), i)) + 1)/2; |
---|
449 | end |
---|
450 | end % X |
---|
451 | |
---|
452 | % create penalized function. Penalize with linear penalty function if |
---|
453 | % violation is severe, otherwise, use exponential penalty. If the |
---|
454 | % 'strict' option has been set, check the constraints, and return INF |
---|
455 | % if any of them are violated. |
---|
456 | function P_fval = funfncP(x, varargin) |
---|
457 | |
---|
458 | % initialize function value |
---|
459 | if (size(x, 2) == 1), P_fval = 0; else P_fval = Nzero.'; end |
---|
460 | |
---|
461 | % initialize x_new array |
---|
462 | x_new = new_x; |
---|
463 | |
---|
464 | % evaluate every vector in x |
---|
465 | for i = 1:size(x, 2) |
---|
466 | |
---|
467 | % reshape x, so it has the same size as the given x0 |
---|
468 | x_new(:) = x(:, i); |
---|
469 | |
---|
470 | % initially, we are optimistic |
---|
471 | linear_eq_Penalty = 0; nonlin_eq_Penalty = 0; |
---|
472 | linear_ineq_Penalty = 0; nonlin_ineq_Penalty = 0; |
---|
473 | |
---|
474 | % Penalize the linear equality constraint violation |
---|
475 | % required: Aeq*x = beq |
---|
476 | if ~isempty(Aeq) && ~isempty(beq) |
---|
477 | lin_eq = abs(Aeq*x_new - beq); |
---|
478 | sumlin_eq = sum(abs(lin_eq(:))); |
---|
479 | if strcmpi(strictness, 'strict') && any(lin_eq > 0) |
---|
480 | P_fval = inf; return |
---|
481 | end |
---|
482 | if sumlin_eq > 200 |
---|
483 | linear_eq_Penalty = exp(200) - 1 + sumlin_eq; |
---|
484 | else |
---|
485 | linear_eq_Penalty = exp(sum(lin_eq)) - 1; |
---|
486 | end |
---|
487 | end |
---|
488 | |
---|
489 | % Penalize the linear inequality constraint violation |
---|
490 | % required: A*x <= b |
---|
491 | if ~isempty(A) && ~isempty(b) |
---|
492 | lin_ineq = A*x_new - b; |
---|
493 | lin_ineq(lin_ineq <= 0) = 0; |
---|
494 | sumlin_ineq = sum(lin_ineq(:)); |
---|
495 | if strcmpi(strictness, 'strict') && any(lin_ineq > 0) |
---|
496 | P_fval = inf; return |
---|
497 | end |
---|
498 | if sumlin_ineq > 200 |
---|
499 | linear_ineq_Penalty = exp(200) - 1 + sumlin_ineq; |
---|
500 | else |
---|
501 | linear_ineq_Penalty = exp(sumlin_ineq) - 1; |
---|
502 | end |
---|
503 | end |
---|
504 | |
---|
505 | % Penalize the non-linear constraint violations |
---|
506 | % required: ceq = 0 and c <= 0 |
---|
507 | if ~isempty(nonlcon) |
---|
508 | [c, ceq] = nonlcon(x_new); %#ok |
---|
509 | if ~isempty(c) |
---|
510 | c = c(:); sumc = sum(c(c > 0)); |
---|
511 | if (strcmpi(strictness, 'strict') || strcmpi(strictness, 'superstrict'))... |
---|
512 | && any(c > ~superstrict*tolCon) |
---|
513 | P_fval = inf; return |
---|
514 | end |
---|
515 | if sumc > 200 |
---|
516 | nonlin_ineq_Penalty = exp200 - 1 + sumc; |
---|
517 | else |
---|
518 | nonlin_ineq_Penalty = exp(sumc) - 1; |
---|
519 | end |
---|
520 | end |
---|
521 | if ~isempty(ceq) |
---|
522 | ceq = abs(ceq(:)); sumceq = sum(ceq(ceq > 0)); |
---|
523 | if (strcmpi(strictness, 'strict') || strcmpi(strictness, 'superstrict'))... |
---|
524 | && any(ceq >= ~superstrict*tolCon) |
---|
525 | P_fval = inf; return |
---|
526 | end |
---|
527 | if sumceq > 200 |
---|
528 | nonlin_eq_Penalty = exp200 - 1 + sumceq; |
---|
529 | else |
---|
530 | nonlin_eq_Penalty = exp(sumceq) - 1; |
---|
531 | end |
---|
532 | end |
---|
533 | end |
---|
534 | |
---|
535 | % evaluate the objective function |
---|
536 | obj_fval = funfcn(x_new, varargin{:}); |
---|
537 | |
---|
538 | % return penalized function value |
---|
539 | P_fval(i) = obj_fval + linear_eq_Penalty + linear_ineq_Penalty + ... |
---|
540 | nonlin_eq_Penalty + nonlin_ineq_Penalty; %#ok |
---|
541 | end |
---|
542 | end % funfncP |
---|
543 | |
---|
544 | % simple wrapper function for output functions; these need to be |
---|
545 | % evaluated with the untransformed variables |
---|
546 | function stop = OutputFcn_wrapper(x, varargin) |
---|
547 | stop = OutputFcn(X(x), varargin{:}); |
---|
548 | end % OutputFcn_wrapper |
---|
549 | |
---|
550 | |
---|
551 | function [output, exitflag] = finalize(x, output, exitflag) |
---|
552 | |
---|
553 | % reshape x so it has the same size as x0 |
---|
554 | x_new = new_x; x_new(:) = x; |
---|
555 | |
---|
556 | % finalize the output |
---|
557 | violated1 = false; violated2 = false; |
---|
558 | violated3 = false; violated4 = false; |
---|
559 | if ~isempty(A) |
---|
560 | Ax = A*x_new; |
---|
561 | violated1 = Ax >= b + tolCon; |
---|
562 | violation = Ax - b; |
---|
563 | violation(~violated1) = 0; |
---|
564 | output.constrviolation.lin_ineq{1} = violated1; |
---|
565 | output.constrviolation.lin_ineq{2} = violation; |
---|
566 | end |
---|
567 | if ~isempty(Aeq) |
---|
568 | Aeqx = Aeq*x_new; |
---|
569 | violated2 = abs(Aeqx - beq) > tolCon; |
---|
570 | violation = Aeqx - beq; |
---|
571 | violation(~violated2) = 0; |
---|
572 | output.constrviolation.lin_eq{1} = violated2; |
---|
573 | output.constrviolation.lin_eq{2} = violation; |
---|
574 | end |
---|
575 | if ~isempty(nonlcon) |
---|
576 | [c, ceq] = nonlcon(x_new); %#ok |
---|
577 | if ~isempty(ceq) |
---|
578 | violated3 = abs(ceq) > tolCon; |
---|
579 | ceq(~violated3) = 0; |
---|
580 | output.constrviolation.nonl_eq{1} = violated3; |
---|
581 | output.constrviolation.nonl_eq{2} = ceq; |
---|
582 | end |
---|
583 | if ~isempty(c) |
---|
584 | violated4 = c > tolCon; c(~violated4) = 0; |
---|
585 | output.constrviolation.nonl_ineq{1} = violated4; |
---|
586 | output.constrviolation.nonl_ineq{2} = c; |
---|
587 | end |
---|
588 | end |
---|
589 | if any([violated1(:); violated2(:); violated3(:); violated4(:)]) |
---|
590 | exitflag = -2; |
---|
591 | output.message = sprintf(... |
---|
592 | ['%s\n Unfortunately, the solution is infeasible for the given value ',... |
---|
593 | 'of OPTIONS.TolCon of %1.6e'], output.message, tolCon); |
---|
594 | else |
---|
595 | if exitflag == 1 |
---|
596 | output.message = sprintf(... |
---|
597 | ['%s\b\n and all constraints are satisfied using ',... |
---|
598 | 'OPTIONS.TolCon of %1.6e.'], output.message, tolCon); |
---|
599 | end |
---|
600 | end |
---|
601 | end % finalize |
---|
602 | end % function |
---|