@Carl Love 

Thanks. Agreed. Perhaps time lapse in posts and I am always confused using the mapleprimes not sure of which post is current or latest.

@Preben Alsholm 

Even in the DAE format one can specify u(0)=alpha, and then add alpha to the parameters list. But alpha should come from the fsolve.

@Preben Alsholm 

If there are multiple solutions for the nonlinear equation, then there will be a different solution for each case. The user will have to decide what makes sense.

If all the roots are meaningul, then there is more than solution for the nonlinear ODE just like we have multiple roots for quadratic or nonlinear equations. 

There is no ODE solver in the world which can solve and give results for all possible solutions numerically (as of now in my experience and search). If the model is known to have multiple solutions, there are approaches to get those solutions and software available for the same (only the final steady state values).

@Carl Love 

The values of the parameters should be the set by the user. Non zero values also work. See code attached. But fsolve needs to be used to find the initial condition.

Download convertimplicitODEtoDAE.mws

@acer 

It works for me in Maple 17 with the commands/code I posted. I was trying to warn the user to not expect this method to work for large scale implicit ODEs.

@Dave L 

I think it will be great to have an fsolve option for a very simple Newton Raphson method, wherein we are guaranteed an ouput with Digits =15  and specified absolute and relative tolerances and a maximum number of iterations. Displaying the iterations is also a great addtion to have so that users will know if there is a progress being made.

@aylin 

Aylin, as shown in the paper, control variables are given by min max commands which are found by solving for the costate equations for bounds if you were to write them.

Once you write down u values, then it is just simulation and is easy to solve. The paper is trying to solve a multi objective problem (by adding weights to different objectives). The problem is badly scaled (in the original Maple code). As of now the major contribution to the objective is from x[1] and not u[1], u[2]. If x[1] is not changing much, just optimizing the other two terms might give faster convergence.

Thanks for this problem, when I find time, I will post a solution for this.

@aylin 

You are welcome. Based on the discussion so far, it is fair to say that you are new to this field.  There are two main approaches to solve these problems, indirect methods and direct methods (Pontyragin). For unconstrained, unbounded problems it is easy to do direct methods. For bounded problems, you have to add costates for each bounds for every state and control variable.

Direct methods can be of 3 types, control vector parameterization, multipleshooting and simultanesous optimizations. If you care only for the solution, then I can write something up and email for the direct methods.

@ 

once we have a meaningful profile for x[1], x[2], x[3] with some optimum constant values for u[1], u[2], one can provide this as an initial guess for the bvp by using approxsoln = sol0 (where sol0 is a good initial guess from the IVP solver).

@aylin 

Is there a physical meaning for x[1], x[2], x[3]? Are they expected to be positive? I wrote a simple code in which I solved only for the original ODEs with u[1] and u[2] as parameters (not changing with time). Increasing u[1] and u[2] seems to increase your objective but x[2] becomes negative.

If x[1], x[2], x[3] are expected to be positive then you need to add adjoint or costate variables for that. You might want to scale the variables and objective so that they are nearly 1. Right now it varies by a wide range.

Download mapleprimeVS.mws

@aylin 

Aylin, thanks for the problem statement and the code. Can you clarify what your objective is (Cost function) and how you got the BCS.

If you are trying to minimize integral of L, then it might result in a singular control problem.

Carl,

Thanks. I asked the original poster to post the original control problem. My current thought is that problem formulation in BVP form is wrong. Typically control variables appear in the original ODEs, not adjoints.

From the Hamiltonian and minimum principles, sometimes one can write control variable as a function of adjoints and state variables (ODE variables).

Depending on the objectives, one of the adjoint variable is likely to be non zero for the boundary condition. Without knowing the original problem statement, we cannot determine this.

To add to this, it is perfectly fine with limitations coming from singularity, etc. It is always left to the user to try to fix these things. 

An important point I would like to bring up is that when CAS is written, please implement the algorithm as is and don't do things that do not reflect the research algorithms or corresponding papers accurately.

http://www.mapleprimes.com/posts/149877-ODEs-PDEs-And-Special-Functions

In the link above, I showed an example of how Maple fails to solve DAEs. The truth of the matter is that at this point in time, Maple does not have a direct DAE solver. It only has ODE solvers. It converts DAE to ODEs and then solves. This point should be made clear in the help page for DAEs. I had spent a lot of time hoping that Maple does what it claims it does. In my opinion, Maplesoft will get more credit if help files reflect and convey these limitations.

All the DAE solvers in Maple have the same issue. The MEBDFI code is based on an algorithm from Dr. Jeff Cash. Unforuntately the Maple code only does ODEs. For DAEs, it converts to ODEs and then solves. This is not the original mebdfi code Professor Cash has published in the literature (implemented in FORTRAN).

My request is that when Maplesoft changes things in algorithms or methods while following a paper or work, please indicate apriori to not waste researchers and students' time. 

PS,  I have modified the Maple code so that one can solve this problem and I can share it if there is interest. I am a longtime user of Maple (from Maple V) and will be a lifetime user of Maple classic worksheet.

Edit:

1. I want to add that Maplesoft folks (in particular Allan Wittkopf) and others are very knowledgeable and very helpful. My guess is they are not staffed enough to update help pages (but this is not an excuse to misrepresent facts). My suggestion would be to add a wiki page for help pages and let the users update them with the caveat that Maplesoft does not own responsbility for the same. This way, when a new user goes to a help page, if there was a linked wiki page then one can read that for typos, algorithms used or modfified.

 2. At this point in time, Maple converts DAEs to explicit ODEs (dy/dt = f). If it is not able to do so, then the code will crash often times sucking up more than 1GB of RAM.

Hi All,

I just installed Maple 18. Yet to check all the new features. Looks like Maple does not provide the option to interact with NAG anymore.

This is a significant backward step. I am not sure why.

Thanks

Venkat Subramanian

@Allan Wittkopf 

Allan,

Thanks for your attention. Just like fsolve/sysnewton, if there are some specific commands that can be used to directly specify equations with variables and initial conditions in proper order (mapped) directly, then that would be great.

Also, please consider providing for direct solution of Mdy/dt =f as discussed here where M is singular for index-1DAE and do direct implicit solver for integration. That will help solve many difficult largescale index-1DAEs.

Note that we should continue to have the current facility to map and identify equations, variables and DAE index,etc. Just that for largescale ( any reasonably sized/difficult) PDEs and DAEs, we need a different/direct approach.

Venkat Subramanian

www.maple.eece.wustl.edu