[4] | 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 |
---|