-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathminimize.m
executable file
·1664 lines (1395 loc) · 65.4 KB
/
minimize.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
%MINIMIZE Solve constrained optimization problems,
% globally or locally
%
% Usage:
% sol = MINIMIZE(func, x0)
% sol = MINIMIZE(..., x0, A, b)
% sol = MINIMIZE(..., b, Aeq, beq)
% sol = MINIMIZE(..., beq, lb, ub)
% sol = MINIMIZE(..., ub, nonlcon)
% sol = MINIMIZE(..., nonlcon, options)
%
% [sol, fval] = MINIMIZE(func, ...)
% [sol, fval, exitflag] = MINIMIZE(func, ...)
% [sol, fval, exitflag, output] = MINIMIZE(func, ...)
% [sol, fval, exitflag, output, grad] = MINIMIZE(func, ...)
% [sol, fval, exitflag, output, gradient, hessian] = MINIMIZE(func, ...)
%
% INPUT ARGUMENTS:
%
% fun, x0 - see FMINSEARCH or FMINLBFGS.
%
% A, b - (OPTIONAL) Linear inequality constraint array and right
% hand side vector.
%
% Aeq, beq - (OPTIONAL) Linear equality constraint array and right
% hand side vector.
%
% lb, ub - (OPTIONAL) lower/upper bound vector or array. Both must have
% the same size as x0.
%
% If no lower bounds exist for one of the variables, then
% supply -inf for that variable. Similarly, if no upper bounds
% exist, supply +inf. If no bounds exist at all, then [lb] and/or
% [ub] may be left empty.
%
% Variables may be fixed in value by setting the corresponding
% lower and upper bounds to exactly the same value.
%
% nonlcon - (OPTIONAL) function handle to general nonlinear constraints,
% inequality and/or equality constraints.
%
% [nonlcon] must return two vectors, [c] and [ceq], containing the
% values for the nonlinear inequality constraints [c] and
% those for the nonlinear equality constraints [ceq] at [x]. (Note:
% these constraints were chosen to be consistent with those of
% fmincon.)
%
% options - (OPTIONAL) an options structure created manually or with
% setoptimoptions().
%
% OUTPUT ARGUMENTS:
%
% sol, fval - the solution vector and the corresponding function value,
% respectively.
%
% exitflag - (See also the help on FMINSEARCH) A flag that specifies the
% reason the algorithm terminated. FMINSEARCH uses only the values
%
% 1 fminsearch converged to a solution x
% 0 Max. # of function evaluations or iterations exceeded
% -1 Algorithm was terminated by the output function.
%
% Since MINIMIZE handles constrained problems, the following
% values were added:
%
% 2 Problem overconstrained by either [lb]/[ub] or
% [Aeq]/[beq] - nothing done
% -2 Problem is infeasible after the optimization (some or
% any of the constraints are violated at the final
% solution).
% -3 INF or NAN encountered during the optimization.
%
% output - (See also the help on FMINSEARCH) A structure that contains
% additional details on the optimization.
%
% gradient/hessian - (only for FMINLBFGS) the gradient and estimated
% hessian at x = sol.
%
% Notes:
%
% If [options] is supplied, then TolX will apply to the transformed
% variables. All other FMINSEARCH parameters should be unaffected.
%
%
% EXAMPLES:
%
% rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2;
%
% >> % Fully unconstrained problem
% >> minimize(rosen, [3 3])
% ans =
% 1.0000 1.0000
%
%
% >> % lower bound constrained
% >> minimize(rosen,[3 3], [],[], [],[], [2 2])
% ans =
% 2.0000 4.0000
%
%
% >> % x(2) fixed at 3
% >> minimize(rosen,[3 3], [],[], [],[], [-inf 3],[inf,3])
% ans =
% 1.7314 3.0000
%
%
% >> % simple linear inequality: x(1) + x(2) <= 1
% >> minimize(rosen,[0; 0], [1 1], 1)
%
% ans =
% 0.6187 0.3813
%
%
% >> % nonlinear inequality: sqrt(x(1)^2 + x(2)^2) <= 1
% >> % nonlinear equality : x(1)^2 + x(2)^3 = 0.5
%
% Execute this m-file:
%
% function test_minimize
% rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2;
%
% options = optimset('TolFun', 1e-8, 'TolX', 1e-8);
%
% minimize(rosen, [3 3], [],[],[],[],[],[],...
% @nonlcon, [], options)
%
% end
% function [c, ceq] = nonlcon(x)
% c = norm(x) - 1;
% ceq = x(1)^2 + x(2)^3 - 0.5;
% end
%
% ans =
% 0.6513 0.4233
%
%
% Of course, any combination of the above constraints is
% also possible.
%
%
% See also: SETOPTIMOPTIONS, FMINSEARCH, FMINLBFGS.
% Please report bugs and inquiries to:
%
% Name : Rody P.S. Oldenhuis
% E-mail : [email protected] (personal)
% [email protected] (professional)
% Affiliation: ispace inc., Japan
% License : BSD 2.0; see license.txt
% FMINSEARCHBND, FMINSEARCHCON and part of the documentation
% for MINIMIZE witten by
%
% Author : John D'Errico
% E-mail : [email protected]
% TODO
%{
WISH: Transformations similar to the ones used for bound constraints should
also be possible for linear constraints; figure this out.
WISH: more checks for inconsistent constraints:
- equality constraints might all lie outside the region defined by the
bounds and linear inequalities
- different linear equality constraints might never intersect
WISH: Include examples and demos in the documentation
WISH: a properly working ConstraintsInObjectiveFunction
FIXME: hess not ouput AT ALL
FIXME: check how functions are evaluated; isn't it better to have objFcn()
and conFcn() do more work?
FIXME: ignore given nonlcon() if ConstraintsInObjectiveFunction is true?
FIXME: fevals is miscounted; possibly fminlbfgs bug
FIXME: TolX, TolFun, TolCon, DiffMinchange, DiffMaxChange should be transformed
FIXME: 'AlwaysHonorConstraints' == 'bounds' does nothing
%}
function [sol, fval,...
exitflag, output,...
grad, hess] = minimize(funfcn, x0,...
A,b,...
Aeq,beq,...
lb,ub,...
nonlcon,...
options, varargin)
% If you find this work useful, please consider a donation:
% https://www.paypal.me/RodyO/3.5
%% Initialization
% Process user input
narg = nargin;
nargo = nargout;
if verLessThan('MATLAB', '8.6')
error(nargchk(2, inf, narg, 'struct')); %#ok<NCHKN>
error(nargchk(0, 6, nargout, 'struct')); %#ok<NCHKM>
else
narginchk (2, inf);
nargoutchk(0, 6);
end
if (narg < 10) || isempty(options), options = setoptimoptions; end
if (narg < 9), nonlcon = ''; end
if (narg < 8), ub = []; end
if (narg < 7), lb = []; end
if (narg < 6), beq = []; end
if (narg < 5), Aeq = []; end
if (narg < 4), b = []; end
if (narg < 3), A = []; end
% Extract options
nonlconFcn_in_objFcn = getoptimoptions('ConstraintsInObjectiveFunction', false);
tolCon = getoptimoptions('TolCon', 1e-8);
algorithm = getoptimoptions('Algorithm', 'fminsearch');
strictness = getoptimoptions('AlwaysHonorConstraints', 'none');
diffMinChange = getoptimoptions('DiffMinChange', 1e-8);
diffMaxChange = getoptimoptions('DiffMaxChange', 1e-1);
finDiffType = getoptimoptions('FinDiffType', 'forward');
OutputFcn = getoptimoptions('OutputFcn', []);
PlotFcn = getoptimoptions('PlotFcn', []);
% Set some logicals for easier reading
have_nonlconFcn = ~isempty(nonlcon) || nonlconFcn_in_objFcn;
have_linineqconFcn = ~isempty(A) && ~isempty(b);
have_lineqconFcn = ~isempty(Aeq) && ~isempty(beq);
create_output = (nargout >= 4);
need_grad = strcmpi(algorithm,'fminlbfgs');
grad_obj_from_objFcn = need_grad && strcmpi(options.GradObj ,'on');
grad_nonlcon_from_nonlconFcn = need_grad && strcmpi(options.GradConstr,'on') && have_nonlconFcn;
do_display = ~isempty(options.Display) && ~strcmpi(options.Display,'off');
do_extended_display = do_display && strcmpi(options.Display,'iter-detailed');
do_global_opt = isempty(x0);
% Do we have an output function?
have_outputFcn = ~isempty(OutputFcn);
have_plotFcn = ~isempty(PlotFcn);
% No x0 given means minimize globally.
if do_global_opt
[sol, fval, exitflag, output, grad, hess] = minimize_globally();
return;
end
% Make copy of UNpenalized function value
UPfval = inf;
% Initialize & define constants
superstrict = false; % initially, don't use superstrict setting
exp50 = exp(50); % maximum penalty
N0 = numel(x0); % variable to check sizes etc.
Nzero = zeros(N0, 1); % often-used zero-matrix
grad = []; % initially, nothing for gradient
sumAT = repmat(sum(A,1).',1,size(b,2)); % column-sum of [A] , transposed and replicated
sumAeqT = repmat(sum(Aeq,1).',1,size(beq,2)); % column-sum of [Aeq], transposed and replicated
% Initialize output structure
output = [];
if create_output
% Fields always present
output.iterations = 0;
output.algorithm = '';
output.message = 'Initializing optimization...';
% Fields present depending on the presence of nonlcon
if ~have_nonlconFcn
output.funcCount = 0;
else
output.ObjfuncCount = 0;
output.ConstrfuncCount = 1; % one evaluation in check_input()
output.constrviolation.nonlin_eq = cell(2,1);
output.constrviolation.nonlin_ineq = cell(2,1);
end
% Fields present depending on the presence of A or Aeq
if have_linineqconFcn
output.constrviolation.lin_ineq = cell(2,1); end
if have_lineqconFcn
output.constrviolation.lin_eq = cell(2,1); end
end
% Save variable with original size
new_x = x0;
% Check for an output/plot functions. If there are any,
% use wrapper functions to call it with un-transformed variable
if have_outputFcn
OutputFcn = options.OutputFcn;
if ~iscell(OutputFcn)
OutputFcn = {OutputFcn}; end
options.OutputFcn = @OutputFcn_wrapper;
end
if have_plotFcn
PlotFcn = options.PlotFcn;
if ~iscell(PlotFcn)
PlotFcn = {PlotFcn}; end
options.PlotFcn = @PlotFcn_wrapper;
end
% Adjust bounds when they are empty
if isempty(lb), lb = -inf(size(x0)); end
if isempty(ub), ub = +inf(size(x0)); end
% Check the user-provided input with nested function check_input
if ~check_input()
return; end
% Force everything to be column vector
ub = ub(:); x0 = x0(:);
lb = lb(:); x0one = ones(size(x0));
% replicate lb or ub when they are scalars, and x0 is not
if isscalar(lb) && (N0 ~= 1), lb = lb*x0one; end
if isscalar(ub) && (N0 ~= 1), ub = ub*x0one; end
% Determine the type of bounds
nf_lb = ~isfinite(lb); nf_ub = ~isfinite(ub);
fix_var = lb == ub; lb_only = ~nf_lb & nf_ub & ~fix_var;
ub_only = nf_lb & ~nf_ub & ~fix_var; unconst = nf_lb & nf_ub & ~fix_var;
lb_ub = ~nf_lb & ~nf_ub & ~fix_var;
%% Optimization
% Force the initial estimate inside the given bounds
x0(x0 < lb) = lb(x0 < lb); x0(x0 > ub) = ub(x0 > ub);
% Transform initial estimate to its unconstrained counterpart
xin = x0; % fixed and unconstrained variables
xin(lb_only) = sqrt(x0(lb_only) - lb(lb_only)); % lower bounds only
xin(ub_only) = sqrt(ub(ub_only) - x0(ub_only)); % upper bounds only
xin(lb_ub) = real(asin( 2*(x0(lb_ub) - lb(lb_ub))./ ...
(ub(lb_ub) - lb(lb_ub)) - 1)); % both upper and lower bounds
xin(fix_var) = [];
% Some more often-used matrices
None = ones(numel(xin)+1,1);
Np1zero = zeros(N0, numel(xin)+1);
% Optimize the problem
try
switch lower(algorithm)
% MATLAB's own derivative-free Nelder-Mead algorithm (FMINSEARCH())
case 'fminsearch'
[presol, fval, exitflag, output_a] = ...
fminsearch(@funfcnT, xin, options);
% Transform solution back to original (bounded) variables...
sol = new_x; sol(:) = X(presol); % with the same size as the original x0
% Evaluate function once more to get unconstrained values
% NOTE: this eval is added to total fevals later
fval(:) = objFcn(sol);
% Steepest descent or Quasi-Newton (limited-memory) BFGS
% (both using gradient) FMINLBFGS(), by Dirk-Jan Kroon
case 'fminlbfgs'
% DEBUG: check gradients
%{
% Numerical
oneHundredth = [
( funfcnT(xin+[1e-2;0])-funfcnT(xin-[1e-2;0]))/2e-2
( funfcnT(xin+[0;1e-2])-funfcnT(xin-[0;1e-2]))/2e-2]
oneTrillionth = [
( funfcnT(xin+[1e-12;0])-funfcnT(xin-[1e-12;0]))/2e-12
( funfcnT(xin+[0;1e-12])-funfcnT(xin-[0;1e-12]))/2e-12]
h = 1e-3;
dfdx = @(f,x) ( (f(x-4*h)-f(x+4*h))/280 + 4*(f(x+3*h)-f(x-3*h))/105 + (f(x-2*h)-f(x+2*h))/5 + 4*(f(x+h)-f(x-h))/5 )/h; %#ok
highOrder = [ dfdx(@(x)funfcnT([x;xin(2)]), xin(1)); dfdx(@(x)funfcnT([xin(1);x]), xin(2)) ]
% As computed in penalized/transformed function
[F,G] = funfcnT(xin)
%}
[presol, fval, exitflag, output_a, ~, hess_a] = ...
fminlbfgs(@funfcnT, xin, options);
% Transform solution back to original (bounded) variables
sol = new_x;
sol(:) = X(presol); % with the same size as the original x0
% Evaluate function some more to get unconstrained values
if grad_obj_from_objFcn
% Function value and gradient
[fval, grad] = objFcn(sol);
if create_output
output_a.funcCount = output_a.funcCount + 1; end
else
% Only function value
[fevals, grad] = computeJacobian(funfcn, sol, objFcn(sol));
if create_output
output_a.funcCount = output_a.funcCount + fevals + 1; end
end
end % switch (algorithm)
catch ME
ME2 = MException([mfilename ':unhandled_error'], [...
'Unhandled problem encountered; please contact ',...
'the author with this exact message, and your ',...
'exact inputs.']);
throw(addCause(ME2,ME));
end
% Copy appropriate fields to the output structure
if create_output
output.message = output_a.message;
output.algorithm = output_a.algorithm;
output.iterations = output_a.iterations;
if ~have_nonlconFcn
output.funcCount = output_a.funcCount + 1;
else
output.ObjfuncCount = output_a.funcCount + 1;
end
end
% Append constraint violations to the output structure, and change the
% exitflag accordingly
[output, exitflag] = finalize(sol, output, exitflag);
% Hessian output
if strcmpi(algorithm, 'fminlbfgs')
% No violation - Hessian estimate equals last-evaluated
if exitflag ~= -2
hess = hess_a;
% Constraints are violated - no unconstrained estimate is available
else
hess = NaN(size(hess_a));
end
end
%% NESTED FUNCTIONS (THE ACTUAL WORK)
% Check user-provided input
function go_on = check_input()
% FIXME: diffMinChange and diffMaxChange can be inconsistent
go_on = true;
% Dimensions & weird input
if (numel(lb) ~= N0 && ~isscalar(lb)) || (numel(ub) ~= N0 && ~isscalar(ub))
error([mfilename ':lb_ub_incompatible_size'], [...
'Size of either [lb] or [ub] incompatible with size ',...
'of [x0].'])
end
if ~isempty(A) && isempty(b)
warning([mfilename 'minimize:Aeq_but_not_beq'], [...
'I received the matrix [A], but you omitted the ',...
'corresponding vector [b].\nI''ll assume a zero-vector ',...
' for [b]...']);
b = zeros(size(A,1), size(x0,2));
end
if ~isempty(Aeq) && isempty(beq)
warning([mfilename ':Aeq_but_not_beq'], [...
'I received the matrix [Aeq], but you omitted the ',...
'corresponding vector [beq].\nI''ll assume a ',...
'zero-vector for [beq]...']);
beq = zeros(size(Aeq,1), size(x0,2));
end
if isempty(Aeq) && ~isempty(beq)
warning([mfilename ':beq_but_not_Aeq'],[...
'I received the vector [beq], but you omitted the ',...
'corresponding matrix [Aeq].\nI''ll ignore the given ',...
'[beq]...']);
beq = [];
end
if isempty(A) && ~isempty(b)
warning([mfilename ':b_but_not_A'], [...
'I received the vector [b], but you omitted the ',...
'corresponding matrix [A].\nI''ll ignore the given ',...
'[b]...']);
b = [];
end
if have_linineqconFcn && size(b,1)~=size(A,1)
error([mfilename ':b_incompatible_with_A'],...
'The size of [b] is incompatible with that of [A].');
end
if have_lineqconFcn && size(beq,1)~=size(Aeq,1)
error([mfilename ':b_incompatible_with_A'],...
'The size of [beq] is incompatible with that of [Aeq].');
end
if ~isvector(x0) && ~isempty(A) && (size(A,2) ~= size(x0,1))
error([mfilename ':A_incompatible_size'], [...
'Linear constraint matrix [A] has incompatible size for ',...
'given [x0].']);
end
if ~isvector(x0) && ~isempty(Aeq) && (size(Aeq,2) ~= size(x0,1))
error([mfilename ':Aeq_incompatible_size'], [...
'Linear constraint matrix [Aeq] has incompatible size ',...
'for given [x0].']);
end
if ~isempty(b) && size(b,2)~=size(x0,2)
error([mfilename ':x0_vector_but_not_b'],[...
'Given linear constraint vector [b] has incompatible ',...
'size with given [x0].']);
end
if ~isempty(beq) && size(beq,2)~=size(x0,2)
error([mfilename ':x0_vector_but_not_beq'], [...
'Given linear constraint vector [beq] has incompatible ',...
'size with given [x0].']);
end
% Functions are not function handles
if ~isa(funfcn, 'function_handle')
error([mfilename ':func_not_a_function'],...
'Objective function must be given as a function handle.');
end
if ~isempty(nonlcon) && ~ischar(nonlcon) && ~isa(nonlcon, 'function_handle')
error([mfilename ':nonlcon_not_a_function'], [...
'Non-linear constraint function must be a function ',...
'handle (advised) or string (discouraged).']);
end
% Check if FMINLBFGS can be executed
if ~isempty(algorithm) && strcmpi(algorithm,'fminlbfgs')
assert(~isempty(which('fminlbfgs')),...
[mfilename ':fminlbfgs_not_present'],[...
'The function FMINLBFGS is not present in the current ',...
'MATLAB path.']);
end
% check output arguments
if nargo >= 5
assert(strcmpi(algorithm,'fminlbfgs'),...
[mfilename ':fminlbfgs_needed_for_derivative_output'], [...
'Gradient or Hessian are only used by the fminlbfgs ',...
'algorithm. Set ''algorithm'' to ''fminlbfgs'' when ',...
'requesting derivative information as an output.']);
end
% evaluate the non-linear constraint function on the initial value,
% to perform initial checks
grad_c = [];
grad_ceq = [];
[c, ceq] = conFcn(x0);
incrementConfuncCount();
% Check sizes of derivatives
if ~isempty(grad_c) && (size(grad_c,2) ~= numel(x0)) && (size(grad_c,1) ~= numel(x0))
error([mfilename ':grad_c_incorrect_size'], [...
'The matrix of gradients of the non-linear in-equality ',...
'constraints\nmust have one of its dimensions equal to ',...
'the number of elements in [x].']);
end
if ~isempty(grad_ceq) && (size(grad_ceq,2) ~= numel(x0)) && (size(grad_ceq,1) ~= numel(x0))
error([mfilename ':grad_ceq_incorrect_size'], [...
'The matrix of gradients of the non-linear equality ',...
'constraints\nmust have one of its dimension equal to ',...
'the number of elements in [x].']);
end
% Test the feasibility of the initial solution (when strict or
% superstrict behavior has been enabled)
if strcmpi(strictness, 'Bounds') || strcmpi(strictness, 'All')
superstrict = strcmpi(strictness, 'All');
if ~isempty(A) && any(any(A*x0 > b))
error([mfilename ':x0_doesnt_satisfy_linear_ineq'],[...
'Initial estimate does not satisfy linear inequality.', ...
'\nPlease provide an initial estimate inside the ',...
'feasible region.']);
end
if ~isempty(Aeq) && any(any(Aeq*x0 ~= beq))
error([mfilename ':x0_doesnt_satisfy_linear_eq'],[...
'Initial estimate does not satisfy linear equality.', ...
'\nPlease provide an initial estimate inside the ',...
'feasible region.']);
end
if have_nonlconFcn
% check [c]
if ~isempty(c) && any(c(:) > ~superstrict*tolCon)
error([mfilename ':x0_doesnt_satisfy_nonlinear_ineq'], [...
'Initial estimate does not satisfy nonlinear ',...
'inequality.\nPlease provide an initial estimate ',...
'inside the feasible region.']);
end
% check [ceq]
if ~isempty(ceq) && any(abs(ceq(:)) >= ~superstrict*tolCon)
error([mfilename ':x0_doesnt_satisfy_nonlinear_eq'], [...
'Initial estimate does not satisfy nonlinear ',...
'equality.\nPlease provide an initial estimate ',...
'inside the feasible region.']);
end
end
end
% Detect and handle degenerate problems
% FIXME: with linear constraints it should be easy to determine if the
% constraints make feasible solutions impossible...
% Impossible constraints
inds = all(A==0,2);
if any(inds) && any(any(b(inds,:) ~= 0))
error([mfilename ':impossible_linear_inequality'],...
'Impossible linear inequality specified.');
end
inds = all(Aeq==0,2);
if any(inds) && any(any(beq(inds,:)~=0))
error([mfilename ':impossible_linear_equality'],...
'Impossible linear equality specified.');
end
% Degenerate constraints
if size(Aeq,1) > N0 && rank(Aeq) == N0
warning([mfilename ':linear_equality_overconstrains'], [...
'Linear equalities overconstrain problem; constrain ',...
'violation is likely.']);
end
if size(Aeq,2) >= N0 && rank(Aeq) >= N0
warning([mfilename ':linear_equality_overconstrains'],...
'Linear equalities define solution - nothing to do.');
sol = Aeq\beq;
fval = objFcn(sol);
new_x = sol;
exitflag = 2;
if create_output
output.iterations = 0;
output.message = 'Linear equalities define solution; nothing to do.';
incrementObjfuncCount();
end
do_display_P = do_display;
do_display = false;
[output, exitflag] = finalize(sol, output, exitflag);
if create_output && exitflag ~= -2
output.message = sprintf(['%s\nFortunately, the solution ',...
'is feasible using OPTIONS.TolCon of %1.6f.'],...
output.message, tolCon);
end
if do_display_P
fprintf(1, output.message); end
go_on = false;
return;
end
% If all variables are fixed, simply return
if sum(lb(:)==ub(:)) == N0
warning([mfilename ':bounds_overconstrain'],...
'Lower and upper bound are equal - nothing to do.');
sol = reshape(lb,size(x0));
fval = objFcn(sol);
new_x = sol;
exitflag = 2;
if create_output
output.iterations = 0;
output.message = 'Lower and upper bound were set equal - nothing to do. ';
incrementObjfuncCount();
end
do_display_P = do_display;
do_display = false;
[output, exitflag] = finalize(sol, output, exitflag);
if create_output && exitflag ~= -2
output.message = sprintf(['%s\nFortunately, the solution ',...
'is feasible using OPTIONS.TolCon of %1.6f.'],...
output.message, tolCon);
end
if do_display_P
fprintf(1, output.message); end
go_on = false;
return;
end
end % check_input
% Counter management
function incrementObjfuncCount(amt)
if create_output
if nargin==0, amt = 1; end
if ~have_nonlconFcn
if ~isfield(output, 'funcCount')
output.funcCount = amt;
else
output.funcCount = output.funcCount + amt;
end
else
if ~isfield(output, 'ObjfuncCount')
output.ObjfuncCount = amt;
else
output.ObjfuncCount = output.ObjfuncCount + amt;
end
end
end
end
function incrementConfuncCount(amt)
if create_output
if nargin==0, amt = 1; end
if ~isfield(output, 'ConstrfuncCount')
output.ConstrfuncCount = amt;
else
output.ConstrfuncCount = output.ConstrfuncCount + amt;
end
end
end
% Evaluate objective function
function varargout = objFcn(x)
[varargout{1:nargout}] = feval(funfcn,...
reshape(x,size(new_x)),...
varargin{:});
end
% Evaluate non-linear constraint function
function [c,ceq, grad_c,grad_ceq] = conFcn(x)
c = [];
ceq = [];
grad_c = [];
grad_ceq = [];
x = reshape(x,size(new_x));
if have_nonlconFcn
if nonlconFcn_in_objFcn
if grad_nonlcon_from_nonlconFcn
if grad_obj_from_objFcn
[~,~, c, ceq, grad_c, grad_ceq] = objFcn(x);
else
[~, c, ceq, grad_c, grad_ceq] = objFcn(x);
end
else
if grad_obj_from_objFcn
[~,~, c, ceq] = objFcn(x);
else
[~, c, ceq] = objFcn(x);
end
end
else
if grad_nonlcon_from_nonlconFcn
[c, ceq, grad_c, grad_ceq] = feval(nonlcon, x, varargin{:});
else
[c, ceq] = feval(nonlcon, x, varargin{:});
end
end
end
end
% Create transformed variable X to conform to upper and lower bounds
function Z = X(x)
% Initialize
if (size(x,2) == 1)
Z = Nzero; rep = 1;
else
Z = Np1zero; rep = None;
end
% First insert fixed values...
y = x0one(:, rep);
y( fix_var,:) = lb(fix_var,rep);
y(~fix_var,:) = x;
x = y;
% ...and transform.
Z(lb_only, :) = lb(lb_only, rep) + x(lb_only, :).^2;
Z(ub_only, :) = ub(ub_only, rep) - x(ub_only, :).^2;
Z(fix_var, :) = lb(fix_var, rep);
Z(unconst, :) = x(unconst, :);
Z(lb_ub, :) = lb(lb_ub, rep) + (ub(lb_ub, rep)-lb(lb_ub, rep)) .* ...
(sin(x(lb_ub, :)) + 1)/2;
end % X
% Derivatives of transformed X
function grad_Z = gradX(x, grad_x)
% ...and compute gradient
grad_Z = grad_x(~fix_var);
grad_Z(lb_only(~fix_var), :) = +2*grad_x(lb_only(~fix_var), :).*x(lb_only(~fix_var), :);
grad_Z(ub_only(~fix_var), :) = -2*grad_x(ub_only(~fix_var), :).*x(ub_only(~fix_var), :);
grad_Z(unconst(~fix_var), :) = grad_x(unconst(~fix_var), :);
grad_Z(lb_ub(~fix_var), :) = grad_x(lb_ub(~fix_var),:).*(ub(lb_ub(~fix_var),:)-lb(lb_ub(~fix_var),:)).*cos(x(lb_ub(~fix_var),:))/2;
end % grad_Z
% Create penalized function. Penalize with linear penalty function if
% violation is severe, otherwise, use exponential penalty. If the
% 'strict' option has been set, check the constraints, and return INF
% if any of them are violated.
function [P_fval, grad_val] = funfcnP(x)
% FIXME: FMINSEARCH() or FMINLBFGS() see this as ONE function
% evaluation. However, multiple evaluations of both objective and nonlinear
% constraint functions may take place
% Initialize function value
if (size(x,2) == 1)
P_fval = 0;
else
P_fval = None.'-1;
end
% Initialize x_new array
x_new = new_x;
% Initialize gradient when needed
if grad_obj_from_objFcn
grad_val = zeros(size(x)); end
% Evaluate every column in x
for ii = 1:size(x,2)
% Reshape x, so it has the same size as the given x0
x_new(:) = x(:,ii);
% Initialize
obj_gradient = 0;
c = []; ceq = [];
grad_c = []; grad_ceq = [];
% Evaluate the objective function, taking care that also
% a gradient and/or contstraint function may be supplied
if grad_obj_from_objFcn
if ~nonlconFcn_in_objFcn
[obj_fval, obj_gradient] = objFcn(x_new);
else
if grad_nonlcon_from_nonlconFcn
arg_out = cell(1, nonlconFcn_in_objFcn+1);
else, arg_out = cell(1, nonlconFcn_in_objFcn+3);
end
[arg_out{:}] = objFcn(x_new);
obj_fval = arg_out{1};
obj_gradient = arg_out{2};
c = arg_out{nonlconFcn_in_objFcn+0};
ceq = arg_out{nonlconFcn_in_objFcn+1};
if grad_nonlcon_from_nonlconFcn
grad_c = arg_out{nonlconFcn_in_objFcn+2};
grad_ceq = arg_out{nonlconFcn_in_objFcn+3};
else
grad_c = ''; % use strings to distinguish them later on
grad_ceq = '';
end
end
objFcn_fevals = 1;
else
% FIXME: grad_nonlcon_from_nonlconFcn?
if ~nonlconFcn_in_objFcn
obj_fval = objFcn(x_new);
else
arg_out = cell(1, nonlconFcn_in_objFcn+1);
[arg_out{:}] = objFcn(x_new);
obj_fval = arg_out{1};
c = arg_out{nonlconFcn_in_objFcn+0};
ceq = arg_out{nonlconFcn_in_objFcn+1};
grad_c = '';
grad_ceq = ''; % use strings to distinguish them later on
end
objFcn_fevals = 1;
if need_grad
[objFevals, obj_gradient] = computeJacobian(funfcn, x_new, obj_fval);
objFcn_fevals = objFcn_fevals + objFevals;
end
end
% Keep track of function evaluations
incrementObjfuncCount(objFcn_fevals);
% Make global copy
UPfval = obj_fval;
% Initially, we are optimistic
linear_eq_Penalty = 0; linear_ineq_Penalty_grad = 0;
linear_ineq_Penalty = 0; linear_eq_Penalty_grad = 0;
nonlin_eq_Penalty = 0; nonlin_eq_Penalty_grad = 0;
nonlin_ineq_Penalty = 0; nonlin_ineq_Penalty_grad = 0;
% Penalize the linear equality constraint violation
% required: Aeq*x = beq
if have_lineqconFcn
lin_eq = Aeq*x_new - beq;
sumlin_eq = sum(abs(lin_eq(:)));
% FIXME: column sum is correct, but does not take into accuont non-violated
% constraints. We'll have to re-compute it
if strcmpi(strictness, 'All') && any(abs(lin_eq) > 0)
P_fval = inf; grad_val = inf; return; end
% FIXME: this really only works with fminsearch; fminlbfgs does not know
% how to handle this...
% compute penalties
linear_eq_Penalty = Penalize(sumlin_eq);
% Also compute derivatives
% (NOTE: since the sum of the ABSOLUTE values is used
% here, the signs are important!)
if grad_obj_from_objFcn && linear_eq_Penalty ~= 0
linear_eq_Penalty_grad = ...
Penalize_grad(sign(lin_eq).*sumAeqT, sumlin_eq);
end
end
% Penalize the linear inequality constraint violation
% required: A*x <= b
if have_linineqconFcn
lin_ineq = A*x_new - b;
lin_ineq(lin_ineq <= 0) = 0;
% FIXME: column sum is correct, but does not take into accuont non-violated
% constraints. We'll have to re-compute it
sumlin_ineq = sum(lin_ineq(:));
sumlin_ineq_grad = sumAT;
if strcmpi(strictness, 'All') && any(lin_ineq > 0)
P_fval = inf; grad_val = inf; return; end
% FIXME: this really only works with fminsearch; fminlbfgs does not know
% how to handle this...
% Compute penalties
linear_ineq_Penalty = Penalize(sumlin_ineq);
% Also compute derivatives
if grad_obj_from_objFcn && linear_ineq_Penalty ~= 0
linear_ineq_Penalty_grad = Penalize_grad(sumlin_ineq_grad, sumlin_ineq); end
end
% Penalize the non-linear constraint violations
% required: ceq = 0 and c <= 0
if have_nonlconFcn
[c, ceq, grad_c, grad_ceq] = conFcn(x_new);
% Central-difference derivatives are computed later;
% the strictness setting might make computing it here
% unneccecary
% Initialize as characters, to distinguish them later on;
% derivatives may be empty, inf, or NaN as returned from [nonlcon]
if isempty(grad_c) , grad_c = ''; end
if isempty(grad_ceq), grad_ceq = ''; end
%{
if ~nonlconFcn_in_objFcn
% Initialize as characters, to distinguish them later on;
% derivatives may be empty, inf, or NaN as returned from [nonlcon]
grad_c = '';
grad_ceq = '';
% Gradients are given explicitly by [nonlcon]
if grad_nonlcon_from_nonlconFcn
[c, ceq, grad_c, grad_ceq] = conFcn(x_new);
% The gradients are not given by [nonlcon]; they have to
% be computed by central differences
else
[c, ceq] = conFcn(x_new);
% Central-difference derivatives are computed later;
% the strictness setting might make computing it here
% unneccecary
end