Advanced Search

Question: Inconsistent behavior with dsolve/numeric and NLPSolve (sqp)

Question:Inconsistent behavior with dsolve/numeric and NLPSolve (sqp)

Dr. Venkat Subramanian 416 Maple 2023
optimization ode dsolve numeric differential-equations nlpsolve
May 31 2023
0

dsolve/numeric + NLPSolve shows inconsistent behavior. This combination is important for parameter estimation and optimal control. Can anyone fix this? Hopefully, I am not making a mistake.

> 

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 dsolve numeric parametric form is buggy and fails for some values of nvar, and works for some values of nvar. Ideally increasing nvar should show convergence with respect to the objective function.

> 

restart:
> 

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):
> 

 

> 

 

> 

ss:=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:
> 

 

> 

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.09025793628011441)

> 

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 = 36893491136404053036)

> 

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

[-.531453381886523, Vector(3, {(1) = .8114345197312305, (2) = 1.189413022736622, (3) = 2.427689509710160})]

> 

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.09025793011810783)

> 

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 = 36893491136437818540)

> 

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

Error, (in Optimization:-NLPSolve) no improved point could be found

> 

 

> 

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.09025792639212991)

> 

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 = 36893491136472713804)

> 

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

[-.537338871783244, Vector(5, {(1) = .7435767224059456, (2) = .9013440906589676, (3) = 1.148772394841535, (4) = 1.629079877401040, (5) = 3.222801229872320})]

> 

 


 

Download test.mw
Please Wait...