@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.
 

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

@vv 

The question was specifically to make NLPSolve work with sqp (gradient approach). 

I am aware of Direct Search package, Global Optimization package from Maple (current version from Optimus and the previous and more robust version which includes the reduced gradient approach (from Pintos) from the early Maple days. 
The code works for me for nvar = 3. I added cmaple to my windows path for a different reason. I am not sure if that made a difference. I am also aware of the analytical solution for a given constant u in a time horizon. 

The goal is to be able to use dsolve numeric/parametric form with sqp approach in Maple. 

The problem statement is given below

dca/dt = -(u+u^2/2)*ca

dcb/dt = u*ca - 0.1*cb

t = 0, ca = 1, cb = 0
Maximize cb(1) subject to the bounds 0 < u< 5.

There is an exact analytical solution for this (Pontyragin's principle). 
This is a reactor problem (with 0.1 replaced by 0) in the second equation introduced/taught by Dr. Larry Biegler at CMU

@vv 

thanks. Matrix form is supposed to be more efficient, but not documented well. 

@tomleslie 

Not sure about your question about nvar,

ss is defined as a procedure which works for any value of nvar = 2,3,4,5 etc

You can ignore the second call for nvar.

If needed you can set nvar:=2 at the top and run the code for just one call. I just showed this so that I can do a loop for nvar and get the objective for different values of nvar as 2,4,8, 16, etc

@minoush82 

I didn't read what happened to your posts, etc.

This is just in support of @acer to make sure he doesn't disappear from this forum.
I often times wondered how he had the time to even help with some very non-significant questions (from my pov), for example, to make a plot prettier, add/modify fonts/shapes, etc.

@acer is technically very good (not just for improving the aesthetics aspects of Maple) and helped me in publishing a peer-reviewed paper for us and is probably one of the very few here who can optimize Maple codes to run at a speed/RAM comparable to C++.
Please check, 
https://www.mapleprimes.com/questions/234050-Will-Anyone-Be-Able-To-Speed-Up-This

FWIW, Maple can compete with even some commercial industrial solvers, I will be posting on the same soon. This wouldn't have been possible without the forum members like @acer, and Maple folks like Allan Wittkopf, Erik Postma, etc

@rlopez 

We recently submitted an arxiv paper for a DAE solver in Maple 

(https://www.mapleprimes.com/questions/235335-Will-Anyone-Be-Able-To-Speed-Up-This)

Arxiv accepts only codes submitted to github. Github currently does not accept Maple files. Perhaps, Maplesoft should move in that direction by interacting with Github?

thanks

The following simple change for finite difference approximation (fourth order as opposed to second order)

for i from 2 to N-1 do
for j from 2 to M-1 do
eq[i,j]:=diff(c[i](t),t)=(16*c[i+1,j](t)-30*c[i,j](t)+16*c[i-1,j](t)-c[i-2,j](t)-c[i+2,j](t))/12/h^2+(16*c[i,j+1](t)-30*c[i,j](t)+16*c[i,j-1](t)-c[i,j-2](t)-c[i,j+2](t))/12/h^2-phi^2*c[i,j](t)^2;
od:od:

helps the code converge faster in x (around N =32 gives 1e-6 error as opposed to N = 64 or 128 for second-order approximation) for example 5. 

Of course, this may not work for all the cases/models (convection, nonlinear diffusivities, singular corners, etc causing unrealistic oscillations), but Maple can be used to develop and try different approximations for spatial derivatives (weak form, collocation, etc).
 

If the solver works or fails for any index-1 DAE that you are testing.