@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

Preben's answer is correct.

Similar problems in cylindrical or other coordinates will have a different particular solution. Typically, separation of variables methods work for non-homogeneous boundary conditions by using c(x,t) = u(x,t) +w(x).

But in this case, that won't work, you can get the answer by using c(x,t) = u(x,t) + w(x) +v(t), and deriving for w(x) and v(t).

More details in the paper, http://www.maple.eece.wustl.edu/pubs/9-JPS-VS-RW-sepvar-2001.pdf. I have also shown this example in the book I wrote on using Maple for chemical engineers.

Good luck. Laplace transform will solve this problem direclty, but you will end up with roots that are repeated and you have to use Heaviside-expansion theorem to invert back to time.

Two flux conditions are common in electrochemical systems (Batteries), and they don't actually have a steady-state. 

Preben's answer is correct.

Similar problems in cylindrical or other coordinates will have a different particular solution. Typically, separation of variables methods work for non-homogeneous boundary conditions by using c(x,t) = u(x,t) +w(x).

But in this case, that won't work, you can get the answer by using c(x,t) = u(x,t) + w(x) +v(t), and deriving for w(x) and v(t).

More details in the paper, http://www.maple.eece.wustl.edu/pubs/9-JPS-VS-RW-sepvar-2001.pdf. I have also shown this example in the book I wrote on using Maple for chemical engineers.

Good luck. Laplace transform will solve this problem direclty, but you will end up with roots that are repeated and you have to use Heaviside-expansion theorem to invert back to time.

Two flux conditions are common in electrochemical systems (Batteries), and they don't actually have a steady-state. 

@ecterrab 

Edgardo, thanks for trying the model I gave. My point is that we should not manipulate index 1 DAE, just take it and solve with implicit formula as explained before.

As you have asked, I am attaching a MATLAB script for the same model (which fails in Maple as of today for N greater than 7). You will see that I am not very good at it. As of now you have to change both N1 and N (both mean the same) for the number of node points. MATLAB solves without any trouble. Plots for both N = 2 and N = 100 are attached. I even ran N = 1000 node points, it ran in 2-3 minutes or so. N = 100 was instantaneous and N = 2 was comparable speed with Maple. In addition, MATLAB file is a pain (I have to print multiple files, merge as a PDF and upload).

For reasons unknown to me MATLAB pushes ode45, but their ode15i and ode15s are their best  solvers in our experience. Note that for this particular model MATLAB is superior, it doesn't mean the same for ill-conditioned problems (Maple's analytical jacobian will make the dsolve/numeric robust and more efficient for ill-conditioned ODEs (small sized) as of today). For index 1 DAEs, I think Maple's way of manipulation and conversion to ODEs should be modified and option to solve directly with M being singular would be great.

MATLABscripts_and_pl.pdfI am sure one can code the same in Mathematica.

Venkat Subramanian

www.maple.eece.wustl.edu

@ecterrab 

Edgardo,

I completely agree with the C code conversion. For example, codegen feature of Maple is a great tool. There is no argument from my side on the successful conversion to C format. My comments are on the implementation of the algorithm, choice of the solver and the method of solving DAEs as ODEs.

When you use standalone C or FORTRAN solvers, very rarely anyone would enter a full dense analytic jacobian for a largescale system. There are more examples that I can provide. For example, Burger's equation in 1 or 2D is a standard problem (solved and demonstrated in Mathematica, MATLAB, etc for example). One will hit the RAM limit as of now with dsolve/numeric when the number of node points is increased.

Please do note that despite the limitation of what we have, our group's first choice is still and will always be Maple to start a model, and when needed hardcoded C or FORTRAN. We very rarely start with other software because for us using Maple helps us develop a very complicated model in 15 minutes to one hour.  However, MATLAB has very good reach among the systems community. For optimization and other toolboxes, MATLAB is far ahead and we use the same for that. In addition, MATLAB allows for files to be sent to non MATLAB users.

Using C directly would take weeks or months to develop these models.

@ecterrab 

Please see a code attached. This comes from finite difference discretization of one PDE with time derivative coupled with one PDE without time derivative. Don't run this code for high values of N (even 5 or 10 would crash any computer), with total number of DAEs less than 25 or so.

The only way to run this code as of today would be to convert all the algebraic equations into ODEs and solve the resulting dy/dt = f system with implicit=true. There again size becomes the limitation.

PDE-DAE-To_post.mw

Hope I am not annoying you, my objective is to continue to restate the current limitation and improve for the future.

Venkat Subramanian

www.maple.eece.wustl.edu

pdsolve/numeric seems to use pdsolve/numeric/sparseLU. I am not sure if is really a sparse solver and if it is available in compiled form.