0001 function t_om_solve_nleqs(quiet)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 if nargin < 1
0013 quiet = 0;
0014 end
0015
0016 core_sp_newton = struct( ...
0017 'alg', 'NEWTON', ...
0018 'name', 'Newton''s', ...
0019 'default_max_it', 10, ...
0020 'need_jac', 1, ...
0021 'update_fcn', @(x, f, J)newton_update_fcn(x, f, J, '') );
0022 core_sp_gs = struct( ...
0023 'alg', 'GS', ...
0024 'name', 'Gauss-Seidel', ...
0025 'default_max_it', 1000, ...
0026 'need_jac', 0, ...
0027 'update_fcn', @(x, f)x_update_fcn2(x, f) );
0028
0029 if have_feature('matlab')
0030
0031 cfg = {
0032 {'DEFAULT', 'default', [] [] },
0033 {'NEWTON', 'Newton', [], [] },
0034 {'FD', 'fast-decoupled Newton',[],[] },
0035 {'FSOLVE', 'fsolve-1', 'fsolve', [] },
0036 {'FSOLVE', 'fsolve-2', 'fsolve', struct('Algorithm', 'trust-region-dogleg') },
0037 {'FSOLVE', 'fsolve-3', 'fsolve', struct('Algorithm', 'trust-region-reflective') },
0038 {'FSOLVE', 'fsolve-4', 'fsolve', struct('Algorithm', 'levenberg-marquardt', 'TolX', 1e-11) },
0039 {'GS', 'Gauss-Seidel',[], [] },
0040 {'CORE-N', 'Newton-CORE', [], core_sp_newton },
0041 {'CORE-GS', 'Gauss-Seidel-CORE',[], core_sp_gs },
0042 };
0043 if have_feature('matlab', 'vnum') <= 7.010
0044 cfg([6]) = [];
0045 end
0046 else
0047
0048 cfg = {
0049 {'DEFAULT', 'default', [] [] },
0050 {'NEWTON', 'Newton', [], [] },
0051 {'FD', 'fast-decoupled Newton',[],[] },
0052 {'FSOLVE', 'fsolve', 'fsolve', struct('TolX', 1e-11) },
0053 {'GS', 'Gauss-Seidel',[], [] },
0054 {'CORE-N', 'Newton-CORE', [], core_sp_newton },
0055 {'CORE-GS', 'Gauss-Seidel-CORE',[], core_sp_gs },
0056 };
0057 end
0058
0059 n = 18;
0060
0061 t_begin(14+n*length(cfg), quiet);
0062
0063 for k = 1:length(cfg)
0064 alg = cfg{k}{1};
0065 name = cfg{k}{2};
0066 check = cfg{k}{3};
0067 opts = cfg{k}{4};
0068 if ~isempty(check) && ~have_feature(check)
0069 t_skip(n, sprintf('%s not installed', name));
0070 else
0071 opt = struct('verbose', 0, 'alg', alg, 'tol', 1e-11);
0072 switch alg
0073 case {'DEFAULT', 'NEWTON'}
0074 case {'FD'}
0075 opt.fd_opt.jac_approx_fcn = @jac_approx_fcn2;
0076 case 'FSOLVE'
0077 opt.fsolve_opt = opts;
0078 case 'GS'
0079 opt.gs_opt.x_update_fcn = @(x, f)x_update_fcn2(x, f);
0080 case {'CORE-N', 'CORE-GS'}
0081 opt.core_sp = opts;
0082 opt.alg = 'CORE';
0083 end
0084
0085 switch alg
0086 case {'DEFAULT', 'NEWTON', 'FSOLVE', 'CORE-N'}
0087 t = sprintf('%s - 2-d function : ', name);
0088 x0 = [-1;0];
0089 om = opt_model;
0090 om.add_var('x', 2, x0);
0091 om.add_nln_constraint('f', 2, 1, @f1, []);
0092 [x, f, e, out, jac] = om.solve(opt);
0093 t_is(e, 1, 12, [t 'success']);
0094 t_is(x, [-3; 4], 8, [t 'x']);
0095 t_is(f, 0, 10, [t 'f']);
0096 switch alg
0097 case {'DEFAULT', 'CORE-N'}
0098 out_alg = 'NEWTON';
0099 otherwise
0100 out_alg = alg;
0101 end
0102 t_ok(strcmp(out.alg, out_alg), [t 'out.alg']);
0103 eJ = [1 1; 6 1];
0104 t_is(jac, eJ, 5.8, [t 'jac']);
0105
0106 t = sprintf('%s - 2-d function (max_it) : ', name);
0107 opt.max_it = 3;
0108 [x, f, e, out] = om.solve(opt);
0109 t_is(e, 0, 12, [t 'no success']);
0110 t_ok(out.iterations == 3 || out.iterations == 4, [t 'iterations']);
0111 opt = rmfield(opt, 'max_it');
0112
0113 t = sprintf('%s - 5-d lin/nonlin function : ', name);
0114 x0_1 = [-1;0];
0115 x0_2 = [0;0;0];
0116 A2 = sparse([2 -1 0; -3 1 -2; 0 5 -4]);
0117 b2 = [-5; 1; -7];
0118 x2 = [-2; 1; 3];
0119 om = opt_model;
0120 om.add_var('x1', 2, x0_1);
0121 om.add_var('x2', 3, x0_2);
0122 om.add_nln_constraint('f', 2, 1, @f1, [], {'x1'});
0123 om.add_lin_constraint('Ax_b', A2, b2, b2, {'x2'});
0124 [x, f, e, out, jac] = om.solve(opt);
0125 t_is(e, 1, 12, [t 'success']);
0126 t_is(x, [-3; 4; -2; 1; 3], 8, [t 'x']);
0127 t_is(f, 0, 10, [t 'f']);
0128 t_ok(strcmp(out.alg, out_alg), [t 'out.alg']);
0129 eJ = [[1 1; 6 1] zeros(2, 3);
0130 zeros(3, 2) A2 ];
0131 t_is(jac, eJ, 5.8, [t 'jac']);
0132 otherwise
0133 t_skip(12, sprintf('not implemented for solver ''%s''', alg));
0134 end
0135
0136 t = sprintf('%s - 2-d function2 (struct) : ', name);
0137 x0 = [1;2];
0138 om = opt_model;
0139 om.add_var('x', 2, x0);
0140 om.add_nln_constraint('f', 2, 1, @f2, []);
0141 [x, f, e] = om.solve(opt);
0142 t_is(e, 1, 12, [t 'success']);
0143 t_is(x, [2; 3], 8, [t 'x']);
0144 t_is(f, 0, 10, [t 'f']);
0145 t_ok(~isfield(om.soln, 'var'), [t 'no parse_soln() outputs']);
0146
0147 opt.max_it = 3;
0148 t = sprintf('%s - 2-d function2 (max_it) : ', name);
0149 [x, f, e, out] = om.solve(opt);
0150 t_is(e, 0, 12, [t 'no success']);
0151 t_ok(out.iterations == 3 || out.iterations == 4, [t 'iterations']);
0152 opt = rmfield(opt, 'max_it');
0153 end
0154 end
0155
0156 t = 'om.soln.';
0157 opt.alg = 'DEFAULT';
0158 x0_1 = [-1;0];
0159 x0_2 = [0;0;0];
0160 A2 = sparse([2 -1 0; -3 1 -2; 0 5 -4]);
0161 b2 = [-5; 1; -7];
0162 x2 = [-2; 1; 3];
0163 om = opt_model;
0164 om.add_var('x1', 2, x0_1);
0165 om.add_var('x2', 3, x0_2);
0166 om.add_nln_constraint('f', 2, 1, @f1, [], {'x1'});
0167 om.add_lin_constraint('Ax_b', A2, b2, b2, {'x2'});
0168 opt.parse_soln = 1;
0169 [x, f, e, out, jac] = om.solve(opt);
0170 t_is(om.soln.x, x, 14, [t 'x']);
0171 t_is(om.soln.f, f, 14, [t 'f']);
0172 t_is(om.soln.eflag, e, 14, [t 'eflag']);
0173 t_ok(strcmp(om.soln.output.alg, out.alg), [t 'output.alg']);
0174 t_is(om.soln.jac, jac, 14, [t 'jac']);
0175
0176 t = 'om.get_soln(''var'', ''x1'') : ';
0177 t_is(om.get_soln('var', 'x1'), x(1:2), 14, [t 'x1']);
0178
0179 t = 'om.get_soln(''var'', ''x'', ''x2'') : ';
0180 t_is(om.get_soln('var', 'x', 'x2'), x(3:5), 14, [t 'x2']);
0181
0182 t = 'om.get_soln(''lin'', ''g'', ''Ax_b'') : ';
0183 g = om.get_soln('lin', 'g', 'Ax_b');
0184 t_is(g{1}, f(3:5), 14, [t 'A * x - u']);
0185 t_is(g{2}, f(3:5), 14, [t 'l - A * x']);
0186
0187 t = 'om.get_soln(''nle'', ''f'') : ';
0188 g = om.get_soln('nle', 'f');
0189 t_is(g, f(1:2), 14, [t 'f']);
0190
0191 t = 'om.get_soln(''nle'', {''g'', ''dg''}, ''f'') : ';
0192 [g, dg] = om.get_soln('nle', {'g', 'dg'}, 'f');
0193 t_is(g, f(1:2), 14, [t 'f']);
0194 t_is(dg, jac(1:2, 1:2), 14, [t 'jac']);
0195
0196 t = 'parse_soln : ';
0197 t_is(om.soln.var.val.x1, om.get_soln('var', 'x1'), 14, [t 'var.val.x1']);
0198 t_is(om.soln.var.val.x2, om.get_soln('var', 'x2'), 14, [t 'var.val.x2']);
0199
0200 t_end;
0201
0202
0203
0204
0205 function [f, J] = f1(x)
0206 if iscell(x)
0207 x = x{1};
0208 end
0209 f = [ x(1) + x(2) - 1;
0210 -x(1)^2 + x(2) + 5 ];
0211 if nargout > 1
0212 J = sparse([1 1; -2*x(1) 1]);
0213 end
0214
0215
0216
0217 function [f, J] = f2(x)
0218 f = [ x(1)^2 + x(1)*x(2) - 10;
0219 x(2) + 3*x(1)*x(2)^2 - 57 ];
0220 if nargout > 1
0221 J = sparse([ 2*x(1)+x(2) x(1);
0222 3*x(2)^2 6*x(1)*x(2)+1 ]);
0223 end
0224
0225 function JJ = jac_approx_fcn2()
0226
0227 J = [7 2; 27 37];
0228 JJ = {J(1,1), J(2,2)};
0229
0230 function x = x_update_fcn2(x, f)
0231
0232 x(1) = sqrt(10 - x(1)*x(2));
0233 x(2) = sqrt((57-x(2))/3/x(1));
0234
0235 function x = newton_update_fcn(x, f, J, lin_solver)
0236 dx = mplinsolve(J, -f, lin_solver);
0237 x = x + dx; 