Venkat Subramanian

386 Reputation

13 Badges

14 years, 235 days

MaplePrimes Activity


These are replies submitted by

Please see the code attached. 
Erik Postma @epostma had helped me in the past to call the MKL solver, so I am just guessing the libraries and trying it. It seems to work. But it may be buggy.

IntelMKL PARDISO can be superior to UMFPACK for very large sparse systems as shown earlier by our group for the paper and codes at https://sites.utexas.edu/maple/phase-field-models-for-li-metal-batteries/


 

NULL

restart:

Digits:=15;

15

(1)

A:=Matrix(4,4,[[-1,2,0,0],[2,-1,2,0],[0,2,-1,2],[0,0,1,-1]],datatype=float[8],storage=sparse);

Matrix(4, 4, {(1, 1) = -1.0, (1, 2) = 2.0, (1, 3) = 0., (1, 4) = 0., (2, 1) = 2.0, (2, 2) = -1.0, (2, 3) = 2.0, (2, 4) = 0., (3, 1) = 0., (3, 2) = 2.0, (3, 3) = -1.0, (3, 4) = 2.0, (4, 1) = 0., (4, 2) = 0., (4, 3) = 1.0, (4, 4) = -1.0})

(2)

b:=Vector(4,[1,0.,0,0],datatype=float[8]):

 

sol:=LinearAlgebra:-LinearSolve(A,b,method=SparseDirectMKL);

Error, invalid input: LinearAlgebra:-LinearSolve expects value for keyword parameter method to be of type identical(none,SparseLU,SparseDirect,SparseDirectMKL,SparseIterative,LU,QR,solve,hybrid,Cholesky,subs,modular), but received SparseDirectMKL

 

s1:=LinearAlgebra:-LUDecomposition(A,method=SparseDirectMKL,output=handle);

before P call
after P call

 

"External/0000029485E225E8"

(3)

s1();

("External/0000029485E225E8")()

(4)


pSolve := define_external("hw_PardisoSolve", ':-MAPLE', ':-LIB' = "linalg"):

x:=Vector(4,datatype=float[8]):

pSolve(b,x,s1);

x;

Vector(4, {(1) = -5.0, (2) = -2.0, (3) = 4.0, (4) = 4.0})

(5)

 

pFactor := define_external("hw_PardisoFactor", ':-MAPLE', ':-LIB' = "linalg"):

pCSF := define_external("hw_SpCompressedSparseForm", ':-MAPLE', ':-LIB' = "linalg"):

 

One can also use pFactor, pCSF, pSolve to get C, R, X and then proceed.

CB,R,X:=pCSF(A,"r",1,1,false,true,false);

CB, R, X := Vector(5, {(1) = 1, (2) = 3, (3) = 6, (4) = 9, (5) = 11}), Vector(10, {(1) = 1, (2) = 2, (3) = 1, (4) = 2, (5) = 3, (6) = 2, (7) = 3, (8) = 4, (9) = 3, (10) = 4}), Vector(10, {(1) = -1., (2) = 2., (3) = 2., (4) = -1., (5) = 2., (6) = 2., (7) = -1., (8) = 2., (9) = 1., (10) = -1.})

(6)

mtype := 11;
handle1 := pFactor(CB, R, X, mtype);

11

 

before P call
after P call

 

"External/0000029485CB0900"

(7)

x:=Vector(4,datatype=float[8]):
pSolve(b,x,handle1):
x;

Vector(4, {(1) = -5.0, (2) = -2.0, (3) = 4.0, (4) = 4.0})

(8)

 

NULL


 

Download buglinearsolve.mw

Such a great thread, thanks all!

@Carl Love 
Not sure I know the difference. I am just looking for something that will call the external integrator without using evalf command.

See a MATLAB call from Maple and the small amount of RAM used (after the MATLAB opens up the first time). (This calls the vdp1 model). The expected wrapper or procedure should call the external integrator like this for a given n(number of odes), f(procedure to find the slopes), initial condition vector (y0), initial time (to) and final time (tf) with default settings for all other parameters

 

@AllanWittkopf is trying to see if he can help. I will update this page. My rough test for a sample code was able to run a model with 50KB RAM consumption for Maple's dsolve/numeric external integrator, this involved guessing the entries for the table creation/modification.

 

restart:

Digits:=15;

15

currentdir();

"C:\Program Files\Maple 2023"

with(Matlab);

[AddTranslator, CloseEngine, ExecuteCommand, FromMFile, FromMatlab, GetVariable, IsAssigned, OpenEngine, SetVariable, chol, closelink, defined, det, dimensions, eig, evalM, fft, getvar, inv, lu, ode15s, ode45, openlink, qr, setvar, size, square, transpose]

ss:=proc(t0,tf,ic0,T,Y)
T,Y := ode45("vdp1", 0 ..tf, ic0);
end proc;

proc (t0, tf, ic0, T, Y) T, Y := Matlab:-ode45("vdp1", 0 .. tf, ic0) end proc

ic0:=[1,0];

[1, 0]

CodeTools:-Usage((ss(0,1,ic0,T,Y)));

memory used=0.79MiB, alloc change=0 bytes, cpu time=0ns, real time=4.29s, gc time=0ns

Vector(57, {(1) = .0, (2) = 0.50237728630191596e-4, (3) = 0.10047545726038319e-3, (4) = 0.1507131858905748e-3, (5) = 0.20095091452076639e-3, (6) = 0.4521395576717244e-3, (7) = 0.7033282008226824e-3, (8) = 0.9545168439736404e-3, (9) = 0.12057054871245985e-2, (10) = 0.24616487028793887e-2, (11) = 0.3717591918634179e-2, (12) = 0.4973535134388969e-2, (13) = 0.622947835014376e-2, (14) = 0.1250919442891771e-1, (15) = 0.1878891050769166e-1, (16) = 0.2506862658646561e-1, (17) = 0.3134834266523956e-1, (18) = 0.5634834266523956e-1, (19) = 0.8134834266523956e-1, (20) = .10634834266523957, (21) = .13134834266523956, (22) = .15634834266523956, (23) = .18134834266523958, (24) = .20634834266523958, (25) = .23134834266523957, (26) = .2563483426652396, (27) = .28134834266523956, (28) = .3063483426652396, (29) = .3313483426652396, (30) = .3563483426652396, (31) = .3813483426652396, (32) = .40634834266523956, (33) = .4313483426652396, (34) = .4563483426652396, (35) = .48134834266523957, (36) = .5063483426652395, (37) = .5313483426652396, (38) = .5563483426652396, (39) = .5813483426652395, (40) = .6063483426652395, (41) = .6313483426652395, (42) = .6563483426652396, (43) = .6813483426652396, (44) = .7063483426652395, (45) = .7313483426652395, (46) = .7563483426652395, (47) = .7813483426652394, (48) = .8063483426652395, (49) = .8313483426652395, (50) = .8563483426652395, (51) = .8813483426652395, (52) = .9063483426652394, (53) = .9313483426652395, (54) = .9485112569989296, (55) = .9656741713326198, (56) = .9828370856663099, (57) = 1.0}, datatype = float[8]), Matrix(57, 2, {(1, 1) = 1.0, (1, 2) = .0, (2, 1) = .9999999987380853, (2, 2) = -0.50237728609061275e-4, (3, 1) = .9999999949523413, (3, 2) = -0.10047545709135335e-3, (4, 1) = .9999999886427678, (4, 2) = -0.15071318532014217e-3, (5, 1) = .9999999798093651, (5, 2) = -0.20095091316873185e-3, (6, 1) = .9999998977849119, (6, 2) = -0.452139542277011e-3, (7, 1) = .9999997526647312, (7, 2) = -0.7033281428979001e-3, (8, 1) = .9999995444488319, (8, 2) = -0.9545166992374083e-3, (9, 1) = .9999992731372271, (9, 2) = -0.12057051955254355e-2, (10, 1) = .999996970144356, (10, 2) = -0.24616462259132695e-2, (11, 1) = .9999930897630864, (11, 2) = -0.37175834032361936e-2, (12, 1) = .9999876319994775, (12, 2) = -0.4973514783111366e-2, (13, 1) = .9999805968620208, (13, 2) = -0.6229438436095478e-2, (14, 1) = .9999217610279643, (14, 2) = -0.12508874315545213e-1, (15, 1) = .9998234934966618, (15, 2) = -0.187878361976329e-1, (16, 1) = .9996857979451632, (16, 2) = -0.25066099711088356e-1, (17, 1) = .9995086794303525, (17, 2) = -0.3134344983295709e-1, (18, 1) = .9984128197702458, (18, 2) = -0.5632104561875044e-1, (19, 1) = .9966928705208215, (19, 2) = -0.8126957756944708e-1, (20, 1) = .9943496689987544, (20, 2) = -.10617984913209057, (21, 1) = .9913842544644557, (21, 2) = -.1310449653737403, (22, 1) = .9877978279054932, (22, 2) = -.1558602995682584, (23, 1) = .9835916754815316, (23, 2) = -.180623437107149, (24, 1) = .9787671017906019, (24, 2) = -.20533409778342732, (25, 1) = .973325392926298, (25, 2) = -.22999409516941732, (26, 1) = .9672677803530982, (26, 2) = -.25460728674271504, (27, 1) = .9605953712287775, (27, 2) = -.27917948363141554, (28, 1) = .9533090911601887, (28, 2) = -.30371835489441507, (29, 1) = .9454096522230412, (29, 2) = -.32823337582573286, (30, 1) = .9368975218313902, (30, 2) = -.35273577439140025, (31, 1) = .9277728623231336, (31, 2) = -.37723843216038444, (32, 1) = .9180354837989904, (32, 2) = -.40175578829756187, (33, 1) = .9076848174098183, (33, 2) = -.4263037870212905, (34, 1) = .896719889317598, (34, 2) = -.4508998266721802, (35, 1) = .8851392702038537, (35, 2) = -.47556266569103106, (36, 1) = .8729410373900864, (36, 2) = -.5003123302542256, (37, 1) = .8601227530027613, (37, 2) = -.5251700638949847, (38, 1) = .8466814429474174, (38, 2) = -.5501582781675165, (39, 1) = .8326135563649153, (39, 2) = -.5753004642541835, (40, 1) = .8179149364435568, (40, 2) = -.6006210932038983, (41, 1) = .8025808031784867, (41, 2) = -.6261455628066843, (42, 1) = .7866057378215437, (42, 2) = -.651900140659666, (43, 1) = .769983652915049, (43, 2) = -.6779118676028639, (44, 1) = .752707773165107, (44, 2) = -.7042084220319103, (45, 1) = .7347706235046282, (45, 2) = -.7308180499230132, (46, 1) = .7161640208762611, (46, 2) = -.7577694811736909, (47, 1) = .6968790574803524, (47, 2) = -.7850917941541237, (48, 1) = .676906097173605, (48, 2) = -.8128141924982113, (49, 1) = .6562347716355753, (49, 2) = -.8409658921908745, (50, 1) = .6348539838487927, (50, 2) = -.8695759788015689, (51, 1) = .6127519109456697, (51, 2) = -.8986731802105471, (52, 1) = .5899160285875202, (52, 2) = -.9282854723793063, (53, 1) = .5663331216220466, (53, 2) = -.9584398809321654, (54, 1) = .5497034312658997, (54, 2) = -.9794686924444002, (55, 1) = .5327104696469083, (55, 2) = -1.0007731387692438, (56, 1) = .5153494429787209, (56, 2) = -1.0223606798647442, (57, 1) = .4976154331426749, (57, 2) = -1.044238264404763}, datatype = float[8])

T[57];Y[57];

1.

Vector[row](%id = 36893490281465971884)

 


 

Download MATLABcall.mw

 

 

Please see the code attached. If anyone can explain the arguments passed then one can possibly write an evalhf-able procedure that will take n (number of odes), f(procedure for the slope), initial time, final time and do all of this with the external integrator and provide the results in evalhf format.
If anyone can show and store the misc and xf called by `dsolve/numeric/SC/IVPrun` this can be done. The goal is to reduce the number of non-evalhf able loops/calls.

(this code is written in Maple 2023).


 

restart:

Digits:=15:

eqs:=diff(ca(t),t)=-(u+u^2/2)*ca(t),diff(cb(t),t)=u*ca(t)-0.1*cb(t);

diff(ca(t), t) = -(u+(1/2)*u^2)*ca(t), diff(cb(t), t) = u*ca(t)-.1*cb(t)

u:=0.1;

.1

infolevel[all]:=10;

10

#soln:=dsolve({eqs,ca(0)=1.0,cb(0)=0.0},type=numeric):

#soln(1);

#showstat(`dsolve/numeric/SC/IVPrun`);

showstat(`dsolve/numeric/SC/hw_int`);


`dsolve/numeric/SC/hw_int` := proc()
local lib;
global `dsolve/numeric/SC/hw_int`;
   1   try
   2       lib := ExternalCalling:-ExternalLibraryName("ode2",'HWFloat');
   3       `dsolve/numeric/SC/hw_int` := define_external('hw_integrate',':-
             MAPLE',(':-LIB') = lib,':-THREAD_SAFE')
       catch :
   4       error
             "'hw_integrate' in external library %1 could not be found/used",
             lib
       end try;
   5   `dsolve/numeric/SC/hw_int`(_passed[1 .. _npassed])
end proc
 

`dsolve/numeric/SC/hw_int` := proc()
local lib;
global `dsolve/numeric/SC/hw_int`;
      try
          lib := ExternalCalling:-ExternalLibraryName("ode2",'HWFloat');
          `dsolve/numeric/SC/hw_int` := define_external('hw_integrate',':-
             MAPLE',(':-LIB') = lib,':-THREAD_SAFE')
       catch :
          error
             "'hw_integrate' in external library %1 could not be found/used",
             lib
       end try;
print(_passed);
      `dsolve/numeric/SC/hw_int`(_passed[1 .. _npassed]);
print("here");
print(_passed);
print(_passed[11]);
end proc;

proc () local lib; global `dsolve/numeric/SC/hw_int`; try lib := ExternalCalling:-ExternalLibraryName("ode2", 'HWFloat'); `dsolve/numeric/SC/hw_int` := define_external('hw_integrate', ':-MAPLE', ':-LIB' = lib, ':-THREAD_SAFE') catch: error "'hw_integrate' in external library %1 could not be found/used", lib end try; print(args); `dsolve/numeric/SC/hw_int`(args[1 .. nargs]); print("here"); print(args); print(args[11]) end proc

soln:=dsolve({eqs,ca(0)=1.0,cb(0)=0.0},type=numeric,output=Array([seq(i/1.*1,i=0..1)])):

dsolve/numeric: entering dsolve/numeric

DEtools/convertsys: Naming is {ca(t) = Y[1] cb(t) = Y[2] diff(ca(t) t) = YP[1] diff(cb(t) t) = YP[2]}

DEtools/convertsys: converted to first-order system Y'(x) = f(x,Y(x))    namely (with Y' represented by YP)

[YP[1] = -.105000000000000*Y[1], YP[2] = .1*Y[1]-.1*Y[2]]

DEtools/convertsys: correspondence between Y[i] names and original functions:

[Y[1] = ca(t), Y[2] = cb(t)]

dsolve/numeric/DAE/rhsproc: Using previously generated rhs proc

dsolve/numeric: the procedure F(x,Y,YP) for computing Y'(x)=f(x,Y(x)) is: proc (N, X, Y, YP) option [Y[1] = ca(t), Y[2] = cb(t)]; YP[1] := -.105000000000000*Y[1]; YP[2] := .1*Y[1]-.1*Y[2]; 0 end proc

dsolve/numeric: initial conditions: x0=0., y0=[1.0, 0.]

dsolve/numeric/SC/preproc: no compilation, _Env_ivpnumeric_ToExternal not set

dsolve/numeric/SC/firststep: Checking ODE for hardware floating point computation

dsolve/numeric/SC/firststep: Initial point OK for hardware floating point computation

dsolve/numeric/SC/firststep: Initial step set to: .0504765875584155

dsolve/numeric/SC/IVPrun: First call, start from initial

dsolve/numeric/SC/IVPrun: Calling external hardware nonstiff integrator with evalhf= true

Vector[row](2, {(1) = 1., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = -.105000000000000, (2) = .100000000000000}), proc (N, X, Y, YP) option `[Y[1] = ca(t), Y[2] = cb(t)]`; YP[1] := (-1)*.105000000000000*Y[1]; YP[2] := .100000000000000*Y[1]+(-1)*.100000000000000*Y[2]; 0 end proc, -1, Vector[row](2, {(1) = .100000000000000, (2) = .100000000000000}), Vector[row](65, {(1) = 2, (2) = 2, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 1, (14) = 2, (15) = 1, (16) = 0, (17) = 0, (18) = 1, (19) = 30000, (20) = 1, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0, (64) = -1, (65) = 0}), Vector[row](28, {(1) = 1., (2) = 1.00000000000000*(1/1000000), (3) = 0.504765875584155e-1, (4) = 5.00001000000000*(1/1000000000000000), (5) = 0., (6) = 0.504765875584155e-1, (7) = 0., (8) = 1.00000000000000*(1/1000000), (9) = 0., (10) = 0., (11) = 0., (12) = 0., (13) = 1., (14) = 0., (15) = .499999999999990, (16) = 0., (17) = 1., (18) = 1., (19) = 0., (20) = 0., (21) = 1., (22) = 1., (23) = 0., (24) = 0., (25) = 1.00000000000000*(1/1000000000000000), (26) = 0., (27) = 0., (28) = 0.}), Array(1..6, 0..2, {(1, 1) = .8758457905848228, (1, 2) = .9121381606660928, (2, 0) = .9121381606660928, (2, 1) = 0.8006442103497206e-1, (2, 2) = 1.0962744218225748, (3, 0) = 1.0962744218225748, (3, 1) = .8912691536438646, (3, 2) = 0.9797584316374096e-1, (4, 0) = 0.9797584316374096e-1, (4, 1) = 1.316703053060327, (4, 2) = .8708776141403536, (5, 0) = .8708776141403536, (5, 1) = .11504702599633022, (5, 2) = 1.5371316842980791, (6, 0) = 1.5371316842980791, (6, 1) = .8509526190302513, (6, 2) = .13130657350916902}, datatype = float[8], order = C_order), 0, Array(1..2, 0..2, {(1, 1) = .0, (1, 2) = 1.0, (2, 0) = 1.0, (2, 1) = .0, (2, 2) = 1.0}, datatype = float[8], order = C_order), 0, 0, 0, 0, 0, 0, 0, 0, Matrix(4, 7, {(1, 1) = 0., (1, 2) = .203125000000000, (1, 3) = .304687500000000, (1, 4) = .750000000000000, (1, 5) = .812500000000000, (1, 6) = .406250000000000, (1, 7) = .812500000000000, (2, 1) = 0.637817382812500e-1, (2, 2) = 0., (2, 3) = .279296875000000, (2, 4) = .272378921508789, (2, 5) = -0.968685150146484e-1, (2, 6) = 0.195693969726562e-1, (2, 7) = .538158416748047, (3, 1) = 0.318908691406250e-1, (3, 2) = 0., (3, 3) = -.343750000000000, (3, 4) = -.335235595703125, (3, 5) = .229614257812500, (3, 6) = .417480468750000, (3, 7) = 11.4807128906250, (4, 1) = 0.971052050590515e-1, (4, 2) = 0., (4, 3) = .403503417968750, (4, 4) = 0.202974677085876e-1, (4, 5) = -0.605428218841553e-2, (4, 6) = -0.477004051208496e-1, (4, 7) = .778585672378540}), Matrix(6, 6, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (1, 6) = 1., (2, 1) = .250000000000000, (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (2, 5) = 0., (2, 6) = 1., (3, 1) = .187500000000000, (3, 2) = .562500000000000, (3, 3) = 0., (3, 4) = 0., (3, 5) = 0., (3, 6) = 2., (4, 1) = .235839843750000, (4, 2) = -.878906250000000, (4, 3) = .890625000000000, (4, 4) = 0., (4, 5) = 0., (4, 6) = .268188476562500, (5, 1) = .127273559570312, (5, 2) = -.500976562500000, (5, 3) = .449218750000000, (5, 4) = -0.128936767578125e-1, (5, 5) = 0., (5, 6) = 0.626220703125000e-1, (6, 1) = -0.927734375000000e-1, (6, 2) = .626220703125000, (6, 3) = -.432617187500000, (6, 4) = .141830444335938, (6, 5) = -0.861053466796875e-1, (6, 6) = .313110351562500}), Vector[row](6, {(1) = 0., (2) = .386000000000000, (3) = .210000000000000, (4) = .630000000000000, (5) = 1., (6) = 1.}), Vector[row](6, {(1) = .250000000000000, (2) = -.104300000000000, (3) = .103500000000000, (4) = -0.362000000000000e-1, (5) = 0., (6) = 0.}), Matrix(6, 5, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (2, 1) = 1.54400000000000, (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (2, 5) = 0., (3, 1) = .946678528081553, (3, 2) = .255701169898258, (3, 3) = 0., (3, 4) = 0., (3, 5) = 0., (4, 1) = 3.31482518706849, (4, 2) = 2.89612401597212, (4, 3) = .998641913997781, (4, 4) = 0., (4, 5) = 0., (5, 1) = 1.22122450922627, (5, 2) = 6.01913448128775, (5, 3) = 12.5370833293209, (5, 4) = -.687886036105895, (5, 5) = 0., (6, 1) = 1.22122450922627, (6, 2) = 6.01913448128775, (6, 3) = 12.5370833293209, (6, 4) = -.687886036105895, (6, 5) = 1.}), Matrix(6, 5, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (2, 1) = -5.66880000000000, (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (2, 5) = 0., (3, 1) = -2.43009335683376, (3, 2) = -.206359915708912, (3, 3) = 0., (3, 4) = 0., (3, 5) = 0., (4, 1) = -.107352905814526, (4, 2) = -9.59456225102190, (4, 3) = -20.4702861480962, (4, 4) = 0., (4, 5) = 0., (5, 1) = 7.49644331396861, (5, 2) = -10.2468043146412, (5, 3) = -33.9999035281991, (5, 4) = 11.7089089320616, (5, 5) = 0., (6, 1) = 8.08324679592241, (6, 2) = -7.98113298806279, (6, 3) = -31.5215943287437, (6, 4) = 16.3193054312314, (6, 5) = -6.05881823883405}), Matrix(3, 5, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (2, 1) = 10.1262350834469, (2, 2) = -7.48799587760763, (2, 3) = -34.8009186155574, (2, 4) = -7.99277170756873, (2, 5) = 1.02513772329562, (3, 1) = -.676280339280690, (3, 2) = 6.08771465167861, (3, 3) = 16.4308432089246, (3, 4) = 24.7672251141837, (3, 5) = -6.59438912571678}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Matrix(2, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0.}), Matrix(2, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Matrix(2, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0.}), Matrix(2, 6, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (1, 6) = 0., (2, 1) = 0., (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (2, 5) = 0., (2, 6) = 0.}), Vector[row](2, {(1) = 0, (2) = 0}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](4, {(1) = 0., (2) = 0., (3) = 0., (4) = 0.}), 0, 0, 0, 0, 0, 0, Vector[row](0, {}), Vector[row](2, {(1) = 0, (2) = 0}), 0

"here"

Vector[row](2, {(1) = .912138160666093, (2) = 0.800644210349721e-1}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = -0.957745068699397e-1, (2) = 0.832073739631121e-1}), proc (N, X, Y, YP) option `[Y[1] = ca(t), Y[2] = cb(t)]`; YP[1] := (-1)*.105000000000000*Y[1]; YP[2] := .100000000000000*Y[1]+(-1)*.100000000000000*Y[2]; 0 end proc, -1, Vector[row](2, {(1) = .100000000000000, (2) = .100000000000000}), Vector[row](65, {(1) = 2, (2) = 2, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 1, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 1, (14) = 2, (15) = 3, (16) = 0, (17) = 0, (18) = 22, (19) = 30000, (20) = 5, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0, (64) = -1, (65) = 0}), Vector[row](28, {(1) = 1., (2) = 1.00000000000000*(1/1000000), (3) = .881714524951008, (4) = 5.00001000000000*(1/1000000000000000), (5) = 0., (6) = 0.504765875584155e-1, (7) = 0., (8) = 1.00000000000000*(1/1000000), (9) = 0., (10) = 0., (11) = 0., (12) = 0., (13) = 1., (14) = 0., (15) = .499999999999990, (16) = 0., (17) = 1., (18) = 1., (19) = 0., (20) = 0., (21) = 1., (22) = 1., (23) = 0., (24) = 0., (25) = 1.00000000000000*(1/1000000000000000), (26) = 0., (27) = 0., (28) = 0.}), Array(1..6, 0..2, {(1, 1) = .8758457905848228, (1, 2) = .9121381606660928, (2, 0) = .9121381606660928, (2, 1) = 0.8006442103497206e-1, (2, 2) = 1.0962744218225748, (3, 0) = 1.0962744218225748, (3, 1) = .8912691536438646, (3, 2) = 0.9797584316374096e-1, (4, 0) = 0.9797584316374096e-1, (4, 1) = 1.316703053060327, (4, 2) = .8708776141403536, (5, 0) = .8708776141403536, (5, 1) = .11504702599633022, (5, 2) = 1.5371316842980791, (6, 0) = 1.5371316842980791, (6, 1) = .8509526190302513, (6, 2) = .13130657350916902}, datatype = float[8], order = C_order), 0, Array(1..2, 0..2, {(1, 1) = .0, (1, 2) = 1.0, (2, 0) = 1.0, (2, 1) = .0, (2, 2) = 1.0}, datatype = float[8], order = C_order), 0, 0, 0, 0, 0, 0, 0, 0, Matrix(4, 7, {(1, 1) = 0., (1, 2) = .203125000000000, (1, 3) = .304687500000000, (1, 4) = .750000000000000, (1, 5) = .812500000000000, (1, 6) = .406250000000000, (1, 7) = .812500000000000, (2, 1) = 0.637817382812500e-1, (2, 2) = 0., (2, 3) = .279296875000000, (2, 4) = .272378921508789, (2, 5) = -0.968685150146484e-1, (2, 6) = 0.195693969726562e-1, (2, 7) = .538158416748047, (3, 1) = 0.318908691406250e-1, (3, 2) = 0., (3, 3) = -.343750000000000, (3, 4) = -.335235595703125, (3, 5) = .229614257812500, (3, 6) = .417480468750000, (3, 7) = 11.4807128906250, (4, 1) = 0.971052050590515e-1, (4, 2) = 0., (4, 3) = .403503417968750, (4, 4) = 0.202974677085876e-1, (4, 5) = -0.605428218841553e-2, (4, 6) = -0.477004051208496e-1, (4, 7) = .778585672378540}), Matrix(6, 6, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (1, 6) = 1., (2, 1) = .250000000000000, (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (2, 5) = 0., (2, 6) = 1., (3, 1) = .187500000000000, (3, 2) = .562500000000000, (3, 3) = 0., (3, 4) = 0., (3, 5) = 0., (3, 6) = 2., (4, 1) = .235839843750000, (4, 2) = -.878906250000000, (4, 3) = .890625000000000, (4, 4) = 0., (4, 5) = 0., (4, 6) = .268188476562500, (5, 1) = .127273559570312, (5, 2) = -.500976562500000, (5, 3) = .449218750000000, (5, 4) = -0.128936767578125e-1, (5, 5) = 0., (5, 6) = 0.626220703125000e-1, (6, 1) = -0.927734375000000e-1, (6, 2) = .626220703125000, (6, 3) = -.432617187500000, (6, 4) = .141830444335938, (6, 5) = -0.861053466796875e-1, (6, 6) = .313110351562500}), Vector[row](6, {(1) = 0., (2) = .386000000000000, (3) = .210000000000000, (4) = .630000000000000, (5) = 1., (6) = 1.}), Vector[row](6, {(1) = .250000000000000, (2) = -.104300000000000, (3) = .103500000000000, (4) = -0.362000000000000e-1, (5) = 0., (6) = 0.}), Matrix(6, 5, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (2, 1) = 1.54400000000000, (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (2, 5) = 0., (3, 1) = .946678528081553, (3, 2) = .255701169898258, (3, 3) = 0., (3, 4) = 0., (3, 5) = 0., (4, 1) = 3.31482518706849, (4, 2) = 2.89612401597212, (4, 3) = .998641913997781, (4, 4) = 0., (4, 5) = 0., (5, 1) = 1.22122450922627, (5, 2) = 6.01913448128775, (5, 3) = 12.5370833293209, (5, 4) = -.687886036105895, (5, 5) = 0., (6, 1) = 1.22122450922627, (6, 2) = 6.01913448128775, (6, 3) = 12.5370833293209, (6, 4) = -.687886036105895, (6, 5) = 1.}), Matrix(6, 5, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (2, 1) = -5.66880000000000, (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (2, 5) = 0., (3, 1) = -2.43009335683376, (3, 2) = -.206359915708912, (3, 3) = 0., (3, 4) = 0., (3, 5) = 0., (4, 1) = -.107352905814526, (4, 2) = -9.59456225102190, (4, 3) = -20.4702861480962, (4, 4) = 0., (4, 5) = 0., (5, 1) = 7.49644331396861, (5, 2) = -10.2468043146412, (5, 3) = -33.9999035281991, (5, 4) = 11.7089089320616, (5, 5) = 0., (6, 1) = 8.08324679592241, (6, 2) = -7.98113298806279, (6, 3) = -31.5215943287437, (6, 4) = 16.3193054312314, (6, 5) = -6.05881823883405}), Matrix(3, 5, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (2, 1) = 10.1262350834469, (2, 2) = -7.48799587760763, (2, 3) = -34.8009186155574, (2, 4) = -7.99277170756873, (2, 5) = 1.02513772329562, (3, 1) = -.676280339280690, (3, 2) = 6.08771465167861, (3, 3) = 16.4308432089246, (3, 4) = 24.7672251141837, (3, 5) = -6.59438912571678}), Vector[row](2, {(1) = .862852234180585, (2) = .121646667791134}), Vector[row](2, {(1) = 0., (2) = 0.}), Vector[row](2, {(1) = -0.873057666308174e-1, (2) = 0.684701246507017e-1}), Vector[row](2, {(1) = 0., (2) = 0.}), Matrix(2, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0.}), Matrix(2, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0.}), Vector[row](2, {(1) = 0., (2) = 0.}), Matrix(2, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0.}), Matrix(2, 6, {(1, 1) = -0.957745068699397e-1, (1, 2) = -0.935578053069745e-1, (1, 3) = -0.925071732935370e-1, (1, 4) = -0.879217450216900e-1, (1, 5) = -0.872861411194737e-1, (1, 6) = -0.914468108408395e-1, (2, 1) = 0.832073739631121e-1, (2, 2) = 0.792621008622733e-1, (2, 3) = 0.774422703214511e-1, (2, 4) = 0.695044462636171e-1, (2, 5) = 0.684014621613813e-1, (2, 6) = 0.755991753809135e-1}), Vector[row](2, {(1) = 0, (2) = 0}), Vector[row](2, {(1) = .912138160666093, (2) = 0.800644210349721e-1}), Vector[row](2, {(1) = .831483491722071, (2) = .146782245215054}), Vector[row](2, {(1) = 1.66780067534233*(1/1000000000), (2) = 6.93290889974829*(1/1000000000)}), Vector[row](2, {(1) = .912138160666093, (2) = 0.800644210349721e-1}), Vector[row](2, {(1) = -0.903900345181292e-1, (2) = 0.737580972195178e-1}), Vector[row](4, {(1) = 0., (2) = 0., (3) = 0., (4) = 0.}), 0, 0, 0, 0, 0, 0, Vector[row](0, {}), Vector[row](2, {(1) = 0, (2) = 0}), 0

Matrix(2, 3, {(1, 1) = 0., (1, 2) = 1.000000000000000, (1, 3) = 0., (2, 1) = 1.000000000000000, (2, 2) = .9003245213933421, (2, 3) = 0.9025791441096123e-1})

dsolve/numeric/SC/IVPrun: Completed external hardware nonstiff integrator

dsolve/numeric/SC/IVPrun: Time for integration= 0.

dsolve/numeric/SC/IVPrun: Time for integration= 0.

 

 

What we need to do is,  write a procedure that will take inputs as the number of odes, n, procedure f, initial conditions for Y vector, initial time, final time, and update the arguments and provide the final Y value from the arguments in evalhf form by direcly calling `dsolve/numeric/SC/hw_int`(_passed[1 .. _npassed]);


 

Download learningdsolvenumeric.mw

Trying to bump up this thread to get the attention of @acer, @tomleslie, and others.
Hopefully, this example convinces everyone that just choosing a different optimization package (global optimization, stochastic, genetic algorithm, etc) may not help in achieving the optimal results if the numerical discretization (in this case nvar for the number time horizons) is not providing sufficient precision. We would need nvar>=256 for this problem and approach.


If anyone has access to old Maple versions, calling the external dsolve/numeric (Without the print statements) and getting an output array as shown in the link below and example towards the end for the NAG subroutine should help.

The link is at 
https://www.maplesoft.com/support/helpJP/Maple/view.aspx?path=NAG/d02cjc

Though Maple does not support NAG toolbox anymore, if anyone has access to Maple 10 or earlier and if you are able to run the example on the NAG help page, please let me know.

Thanks

 

 

@sursumCorda 

This is based on my limited testing of combining Maple's code with its global optimization toolbox.

Currently, you pay a separate license for the same.

Currently, this is offered by Nimbus which provides Genetic Algorithms.
A decade or so earlier, Maple provided the Global Optimization toolbox from Pintur. That had method=reducedgradient which was very robust and fast. 
Hope this helps. 

@tomleslie 

I like DirectSearch and I have nothing against it or genetic algorithms (GA) which are heavily used in the literature.

This discussion and analysis (from @acer and the solution from him) is the best and most efficient way to solve this model. 


The next biggest step can be wrapping dsolve/numeric (with or without compile=true) and making the NLPSolve work with the same in evalhf form for computational efficiency. 

The last 5-10 years have seen an exponential increase in open access/open software codes in Python, Julia, etc. Many times I have wondered why students/researchers used GA for well-defined models (with convexity). Now, I partially understand the reason. Many times wrappers are written for optimizers or different algorithms and then called from user-friendly software. But the wrapper is probably poorly written so that the gradient information is not properly passed with a sufficient number of options (tolerances in this case).

PS, note that only gradient-based methods have theoretical backing for locally optimal results for optimization and I am weary of using methods that don't have a good theoretical basis, though they might work very well (but they will be always slower for problems that can be solved with gradient based methods).

 

@acer 

I think any opimizer or nonlinearsolver should take the required tolerance as input and return the current-best objective, iterate.
@acer, I can confirm that this approach works with nonlinear constraints + dsolve/numeric/parameters + sqp. thanks again

@acer 

There is no point in doing optimization to 1e-16 accuracy when we are trying to increase nvar(number of stages), solving ODEs/DAEs to 1e-6 accuracy.

I like your answer and it works.

Note that for some weird reason, for this problem stiff=true in the call to dsolve/numeric works without tweaking anything in the dsolve command or optimization command.  

Knowing where it fails can help. For example, ideally NLPSolve should have given the answer that "objective was solved to 1e-6 precision, but not beyond that", instead of just returning "initial point satisfies first-order optimality condition".

fsolve has the same issue as well, and if I am not wrong, we can't force it to solve to a prescribed precision.

@acer, as again, thanks lot. I believe that you can fix almost any issues!

Adding abstol=1e-12, reltol=1e-12 to dsolve/numeric calls solves the problem for this model. (Thanks to Allan Wittkopf)

He had also suggested using Digits:=18 inside ss (I am not a big fan of this)

With default tolerances use of, maxstep = 0.05 (assuming that your maximum time horizon is 0.1 with nvar = 10)seems to work as well.

fdiff or the procedures used by NLPSolve should be made more robust.

@dharr 

Your code is not using dsolve numeric in parametric form. I will be posting the use of parametric form working with NLPSolve (thanks to Allan Wittkopf).

As I said, we don't want analytical solutions.

Also, one parameter optimization is different from multiple parameters (scalar vs vector).

To be clear, I am not looking for an answer for the model. My goal was to use NLPSolve with dsolve/numeric/parameters with sqp. 
 

@dharr 

thanks

When dsolve numeric is replaced with the RK2 method, the code works. Typically, dsolve numeric will be faster than writing for loops to do numerical integration of ODEs.

This is why I think that the issue is with dsolve/numeric.
 

restart:

 

Test code written by Dr. Venkat Subramanian at UT Austin, 05/31/2023. This code uses CVP approach (piecwise constant) to perform optimal control. NLPSolve combination with RK2 approach to integrate ODEs (constant step size) works well for this problem.

Digits:=15;

15

 

eqodes:=[diff(ca(t),t)=-(u+u^2/2)*1.0*ca(t),diff(cb(t),t)=1.0*u*ca(t)-0.1*cb(t)];

[diff(ca(t), t) = -1.0*(u+(1/2)*u^2)*ca(t), diff(cb(t), t) = 1.0*u*ca(t)-.1*cb(t)]

soln:=dsolve({op(eqodes),ca(0)=alpha,cb(0)=beta},type=numeric,'parameters'=[alpha,beta,u],compile=true,savebinary=true):

RK2:=proc(alpha,beta,NN,u,dt)#this procedure can be made efficient in vector form
local j,c1mid,c2mid,c10,c20,c1,c2;
c10:=alpha;c20:=beta;
for j from 1 to NN do
  c1mid:=c10+dt/NN*(-(u+u^2/2)*c10):
  c2mid:=c20+dt/NN*(u*c10)-0.1*dt/NN*c20:
  c1:=c10/2+c1mid/2+dt/NN/2*(-(u+u^2/2)*c1mid):
  c2:=c20/2+c2mid/2+dt/NN/2*(u*c1mid)-0.1*dt/NN/2*c2mid:
  c10:=c1:c20:=c2:
  od:
[c10,c20];
end proc:

soln('parameters'=[1,0,0.1]);soln(0.1);

[alpha = 1., beta = 0., u = .1]

[t = .1, ca(t) = HFloat(0.9895549324188543), cb(t) = HFloat(0.009898024276129616)]

RK2(1,0,20,0.1,0.1);#reasonable accuracy with just 20 time steps of RK2

[.989554933045397, 0.989802232901406e-2]

ss0:=proc(x)
interface(warnlevel=0):
#if  type(x[1],numeric)
if  type(x,Vector)
then

local z1,n1,i,c10,c20,dt,u;
global soln,nvar;
dt:=evalf(1.0/nvar):
c10:=1.0:c20:=0.0:
for i from 1 to nvar do
u:=x[i]:
soln('parameters'=[c10,c20,u]):
z1:=soln(dt):
c10:=subs(z1,ca(t)):c20:=subs(z1,cb(t)):
od:
-c20;
 else 'procname'(args):

end if:

end proc:

ss:=proc(x)#based on RK2
interface(warnlevel=0):
#if  type(x[1],numeric)
if  type(x,Vector)
then

local z1,n1,i,c10,c20,dt,u,NN;
global soln,nvar;
dt:=evalf(1.0/nvar):
NN:=20;
c10:=1.0:c20:=0.0:
for i from 1 to nvar do
u:=x[i]:
z1:=RK2(c10,c20,NN,u,dt):
c10:=z1[1]:c20:=z1[2]:
od:
-c20;
 else 'procname'(args):

end if:

end proc:

nvar:=3;#code works for nvar:=3, but not for nvar:=2

3

ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float[8]);

Vector(3, {(1) = .1000000000000000, (2) = .1000000000000000, (3) = .1000000000000000})

 

ss(ic0);

HFloat(-0.09025778155141699)

bl := Vector(nvar,[seq(0.,i=1..nvar)]);bu := Vector(nvar,[seq(5.,i=1..nvar)]);

Vector(3, {(1) = 0., (2) = 0., (3) = 0.})

Vector[column](%id = 36893490450763924652)

infolevel[all]:=1;

1

Optimization:-NLPSolve(nvar,ss,[],initialpoint=ic0,[bl,bu]);

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 3

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

[-.531367180853919, Vector(3, {(1) = .8116892515974860, (2) = 1.189817750741097, (3) = 2.426110609252688})]

nvar:=2;#code works for nvar:=3, but not for nvar:=2

2

ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float[8]);

Vector(2, {(1) = .1000000000000000, (2) = .1000000000000000})

 

ss(ic0);

HFloat(-0.09025762200103044)

bl := Vector(nvar,[seq(0.,i=1..nvar)]);bu := Vector(nvar,[seq(5.,i=1..nvar)]);

Vector(2, {(1) = 0., (2) = 0.})

Vector[column](%id = 36893490450763939460)

infolevel[all]:=1;

1

Optimization:-NLPSolve(nvar,ss,[],initialpoint=ic0,[bl,bu]);

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 2

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

[-.523759339206159, Vector(2, {(1) = .9103276760025985, (2) = 1.925973211795306})]

 

nvar:=5;#code works for nvar:=3, but not for nvar:=2

5

ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float[8]);

Vector(5, {(1) = .1000000000000000, (2) = .1000000000000000, (3) = .1000000000000000, (4) = .1000000000000000, (5) = .1000000000000000})

 

ss(ic0);

HFloat(-0.09025786314622318)

bl := Vector(nvar,[seq(0.,i=1..nvar)]);bu := Vector(nvar,[seq(5.,i=1..nvar)]);

Vector(5, {(1) = 0., (2) = 0., (3) = 0., (4) = 0., (5) = 0.})

Vector[column](%id = 36893490450873412892)

infolevel[all]:=1;

1

Optimization:-NLPSolve(nvar,ss,[],initialpoint=ic0,[bl,bu]);

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 5

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

[-.537292701143927, Vector(5, {(1) = .7437459723532210, (2) = .9015469842007691, (3) = 1.149056357665623, (4) = 1.629519387333676, (5) = 3.220718773451731})]

infolevel[all]:=0:

#Note that results from nvar =2 are good initial guesses for nvar =4, etc, but it is not implemented in this code.

# objective is approaching convergence with 16 stages for control variable.

for j from 1 to 4 do
nvar:=2^j:
ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float[8]);
bl := Vector(nvar,[seq(0.,i=1..nvar)]);bu := Vector(nvar,[seq(5.,i=1..nvar)]);
obj[j]:=Optimization:-NLPSolve(nvar,ss,[],initialpoint=ic0,[bl,bu]):
print(2^j,obj[j][1]):
gc();
od:

2, -.523759339206158692

4, -.535093353250732817

8, -.540523213471601371

16, -.542691160681018525

[seq(obj[4][2][i],i=1..16)];#Control values

[HFloat(0.6844327854667692), HFloat(0.7213390420309091), HFloat(0.762134925769591), HFloat(0.8075994009309487), HFloat(0.8587549491734409), HFloat(0.9169717511913674), HFloat(0.9841099356697114), HFloat(1.0627974101121582), HFloat(1.1569237891930957), HFloat(1.2724091321519309), HFloat(1.4190145332788868), HFloat(1.6140310681691525), HFloat(1.89202739512972), HFloat(2.3355791576022167), HFloat(3.214921631547972), HFloat(5.0)]

 


 

Download testRk2.mw

@dharr 

I am a long-time Maple user, and I would appreciate it if you can provide any code/answer without the shortcuts/etc as much as possible. Also, for any procedures please don't add any inputs/variables/descriptions/variable types/ etc unless they are critically needed to make the code work for that purpose.

One of the charms of Maple is its symbolic capability, and keeping the codes/commands minimal and getting the results we seek are the best options for people like me who still miss the classic worksheet and old Maple GUIs.

Thanks

Venkat

 

@dharr 

As stated early, I am only interested in the sqp and gradient-based method. Nonlinear simplex cannot handle bounds.
My interest lies in sqp as I have to be able to able to change nvar as input to a procedure to predict the objective, and possibly nonlinear constraints as an objective

1 2 3 4 5 6 7 Last Page 1 of 17