DAE daspk solver problem

classic Classic list List threaded Threaded
15 messages Options
Reply | Threaded
Open this post in threaded view
|

DAE daspk solver problem

Lazy_Tom
Hi all,

I'd like to use the daspk solver for DAE that comes with Octave but I can't find any official example on how to use it so I decided to give it a go myself.
I've tried to solve the pendulum equations in cartesian coordinates but I'm having a hard time working out why the solver it is not working.

the variables are as follows

x(1)=position x
x(2)=velocity x
x(3)=position y
x(4)=velocity y
x(5)=reaction force intensity

xdot(1)=velocity x
xdot(2)=acceleration x
xdot(4)=velocity y
xdot(5)=velocity y

and the script is

t=linspace(0,5,1000)';
m=1;
g=9.81;

function f=effe(x,xdot, t)

m=1;
g=9.81;

f(1) = xdot(1)-x(2);
f(2) = m*xdot(2)+x(5)*x(1);
f(3) = xdot(3)-x(4);
f(4) = m*xdot(4) + x(5)*x(3) +m*g;
f(5) = x(1)^2 + x(3)^2 - 1;
endfunction
daspk_options('print initial condition info', 1);
daspk_options('compute consistent initial condition', 1);
daspk_options('algebraic variables', [0,0,0,0,1]);

[x, XDOT, ISTATE, MSG] = daspk (@effe, [1,0, 0, 0, m*g], [0,0,0,g,0], t);

any help would be appreciated.

Regards

Marco
Reply | Threaded
Open this post in threaded view
|

Re: DAE daspk solver problem

Sebastian Schöps
Lazy_Tom wrote
f(1) = xdot(1)-x(2);
f(2) = m*xdot(2)+x(5)*x(1);
f(3) = xdot(3)-x(4);
f(4) = m*xdot(4) + x(5)*x(3) +m*g;
f(5) = x(1)^2 + x(3)^2 - 1;
I think there are a few typos, e.g., it should be -x(5)*x(3) in eq (4) but this is probably not the issue here. You are trying to solve a very difficult problem. This formulation of the pendulum is a DAE of index 3 and afaik daspk can only deal with DAEs of index<3. There are many good books on this topic. If you are looking for a good free ressource, you may want to look at https://www.mathematik.hu-berlin.de/~steffen/pub/introduction_to_daes_497.pdf (page 2 and 13).

Seb.
Reply | Threaded
Open this post in threaded view
|

Re: DAE daspk solver problem

Lazy_Tom
Sebastian Schöps wrote
Lazy_Tom wrote
f(1) = xdot(1)-x(2);
f(2) = m*xdot(2)+x(5)*x(1);
f(3) = xdot(3)-x(4);
f(4) = m*xdot(4) + x(5)*x(3) +m*g;
f(5) = x(1)^2 + x(3)^2 - 1;
I think there are a few typos, e.g., it should be -x(5)*x(3) in eq (4) but this is probably not the issue here. You are trying to solve a very difficult problem. This formulation of the pendulum is a DAE of index 3 and afaik daspk can only deal with DAEs of index<3. There are many good books on this topic. If you are looking for a good free ressource, you may want to look at https://www.mathematik.hu-berlin.de/~steffen/pub/introduction_to_daes_497.pdf (page 2 and 13).

Seb.
Hi Seb,

thanks for pointing out daspk max index, I wasn't aware. it spared a lot of headache.

indeed, for some obscure reason I was convinced that the option "maximum order" was actually "maximum index"... -.-

I don't know why DAE solvers index capabilities are not reported in Octave manual. Matlab at least tells you that you can go no further than index 1 (at least without doing some "magic"...)

about the paper you cited, yup, I'm aware of it and I agree, it is quite nice.

Well, at least I'll gave it a try, I think that for DAE problems I'll move to openModelica

Cheers,

Marco


Reply | Threaded
Open this post in threaded view
|

Re: DAE daspk solver problem

Sebastian Schöps
Lazy_Tom wrote
I don't know why DAE solvers index capabilities are not reported in Octave manual. Matlab at least tells you that you can go no further than index 1 (at least without doing some "magic"...)
You are right. Feel free to oben a bug report at savannah, maybe immediately suggesting a string that can be added to the help. However, I assume that dassl and daspk will be deprecated in some future version since we are moving to odeXY compatible solvers.

Lazy_Tom wrote
about the paper you cited, yup, I'm aware of it and I agree, it is quite nice.
Well, at least I'll gave it a try, I think that for DAE problems I'll move to openModelica
I think that ode5r from odepkg should be able to solve index3 problems since its based on Radau 5 [1]. We hope to have a new release of odepkg quickly after Octave 4.2. I will keep in mind to include such information in the help.

Regarding openModelica: yes, they do some clever tricks for higher index. As far as I recall it's based on dummy derivatives [2]. Such stuff is easier for modelica since it is in some sense a symbolic language. However, to my best knowledge, it is not really well suited for big problems. So, as usual, the best choice depends on what you want to do :)

Seb.

[1] http://www.unige.ch/~hairer/prog/stiff/radau5.f
[2] http://dx.doi.org/10.1137/0914043
Reply | Threaded
Open this post in threaded view
|

Re: DAE daspk solver problem

Bill Greene-3
I was able to solve this using ode5r from odepkg. My test code is below:

function radau5PendulumDAE

vopt = odeset ('InitialStep', 1e-3, 'AbsTol', 1e-8);
len=1; theta0=pi/2; m=1; g=9.81;
y0=[ len*sin(theta0), -len*cos(theta0), 0, 0, m*g];
odef = @(t,x) odeFunc (t, x, m, g, len)
sol=ode5r (odef, [0, 5], y0, vopt);
figure; plot(sol.x, sol.y(:,1));
end

function y = odeFunc (t, x, m, g, len)
lam = x(5);
y=[x(2); -lam*x(1); x(4); -lam*x(3)-m*g; x(1)^2 + x(3)^2 - 1];
end

On Fri, Oct 28, 2016 at 2:18 AM, Sebastian Schöps <[hidden email]> wrote:
Lazy_Tom wrote
> I don't know why DAE solvers index capabilities are not reported in Octave
> manual. Matlab at least tells you that you can go no further than index 1
> (at least without doing some "magic"...)

You are right. Feel free to oben a bug report at savannah, maybe immediately
suggesting a string that can be added to the help. However, I assume that
dassl and daspk will be deprecated in some future version since we are
moving to odeXY compatible solvers.


Lazy_Tom wrote
> about the paper you cited, yup, I'm aware of it and I agree, it is quite
> nice.
> Well, at least I'll gave it a try, I think that for DAE problems I'll move
> to openModelica

I think that ode5r from odepkg should be able to solve index3 problems since
its based on Radau 5 [1]. We hope to have a new release of odepkg quickly
after Octave 4.2. I will keep in mind to include such information in the
help.

Regarding openModelica: yes, they do some clever tricks for higher index. As
far as I recall it's based on dummy derivatives [2]. Such stuff is easier
for modelica since it is in some sense a symbolic language. However, to my
best knowledge, it is not really well suited for big problems. So, as usual,
the best choice depends on what you want to do :)

Seb.

[1] http://www.unige.ch/~hairer/prog/stiff/radau5.f
[2] http://dx.doi.org/10.1137/0914043



--
View this message in context: http://octave.1599824.n4.nabble.com/DAE-daspk-solver-problem-tp4680377p4680395.html
Sent from the Octave - General mailing list archive at Nabble.com.

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: DAE daspk solver problem

sshah
In reply to this post by Sebastian Schöps
As I understand it, Radau5.f demands a structure of DAEs with Linearly implicit form:

M Y' = F(Y,T)

DASPK admits allows a more general form F(Y,Y',T)=0.  There can be problems that can not be transformed to the linearly implicit form.  

For large scale problems, (> 1e5 variables), retaining sparsity is an important issue, so DASPK structure can offer some advantages.  Even when the equations can be symbolically transformed, symbolic manipulation may not scale well (perhaps the problem with scaling in Modelica),  or may destroy sparsity and therefore have inefficient discrete propagation.

Lastly, the fortran DASPK / DASKR can take advantage of implicit sparse solvers.  Do octave implementations allow this?



Reply | Threaded
Open this post in threaded view
|

Re: DAE daspk solver problem

Bill Greene-3
I did not try this simple test in daspk but I did try it in ida (with C-api).
Even after carefully
insuring that I had consistent initial conditions, ida could not advance to the
first time step. So I am very doubtful that daspk will solve this example in its
index-3 form.

Frankly, I was quite surprised that radau5 was able to do so particularly since
it doesn't appear to have an option to specify an initial y'. 


On Fri, Oct 28, 2016 at 1:17 PM, SunilShah <[hidden email]> wrote:
As I understand it, Radau5.f demands a structure of DAEs with Linearly
implicit form:

M Y' = F(Y,T)

DASPK admits allows a more general form F(Y,Y',T)=0.  There can be problems
that can not be transformed to the linearly implicit form.

For large scale problems, (> 1e5 variables), retaining sparsity is an
important issue, so DASPK structure can offer some advantages.  Even when
the equations can be symbolically transformed, symbolic manipulation may not
scale well (perhaps the problem with scaling in Modelica),  or may destroy
sparsity and therefore have inefficient discrete propagation.

Lastly, the fortran DASPK / DASKR can take advantage of implicit sparse
solvers.  Do octave implementations allow this?







--
View this message in context: http://octave.1599824.n4.nabble.com/DAE-daspk-solver-problem-tp4680377p4680409.html
Sent from the Octave - General mailing list archive at Nabble.com.

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: DAE daspk solver problem

sshah
DASPK/DASKR are not expected to solve index 2+ problems like the index-3 pendulum problem.

My comment is about problem structure (mass matrix vs fully implicit), sparsity and scaling.
Reply | Threaded
Open this post in threaded view
|

Re: DAE daspk solver problem

kingcrimson
In reply to this post by sshah

On 28 Oct 2016, at 19:17, SunilShah <[hidden email]> wrote:

> For large scale problems, (> 1e5 variables), retaining sparsity is an
> important issue, so DASPK structure can offer some advantages.  Even when
> the equations can be symbolically transformed, symbolic manipulation may not
> scale well (perhaps the problem with scaling in Modelica),  or may destroy
> sparsity and therefore have inefficient discrete propagation.

Unfortunately the sparsity is not exploited in the Octave interface.

> Lastly, the fortran DASPK / DASKR can take advantage of implicit sparse
> solvers.  Do octave implementations allow this?

No, AFAIK there is no option in Octave to invoke the matrix-free iterative solver
in DASPK from Octave.

c.


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: DAE daspk solver problem

Lazy_Tom
In reply to this post by Bill Greene-3
Bill Greene-3 wrote
I was able to solve this using ode5r from odepkg. My test code is below:

function radau5PendulumDAE

vopt = odeset ('InitialStep', 1e-3, 'AbsTol', 1e-8);
len=1; theta0=pi/2; m=1; g=9.81;
y0=[ len*sin(theta0), -len*cos(theta0), 0, 0, m*g];
odef = @(t,x) odeFunc (t, x, m, g, len)
sol=ode5r (odef, [0, 5], y0, vopt);
figure; plot(sol.x, sol.y(:,1));
end

function y = odeFunc (t, x, m, g, len)
lam = x(5);
y=[x(2); -lam*x(1); x(4); -lam*x(3)-m*g; x(1)^2 + x(3)^2 - 1];
end
Hi Bill,

I don't think your description is correct, at the moment you are solving a regular ODE. The solver is not aware that the last equation doesn't contain the derivative of x(5). And this needs to be specified with odeset 'Mass' option and other parameters

Out of curiosity, getting odest documentation (command 'doc odepkg' doesn't work on my windows release) and odepkg.pdf document has been a problem. The pdf is mentioned in odeset help, odepkg help and website (http://octave.sourceforge.net/odepkg/overview.html) but to the best of my knowledge the actual location of it doesn't exist. In fact, after installing the odepkg I did a search odepkg.pdf on my octave directory with no results. What I've found though is an odepkg.texi file which is indeed the package manual but it's not a pdf.

This is how I finally found out about the 'Mass' option (yes there is a tiny example in odepkg_examples_dae but I don't think this is enough)

Do you or anyone else have this problem as well?

Marco
Reply | Threaded
Open this post in threaded view
|

Re: DAE daspk solver problem

Sebastian Schöps
Lazy_Tom wrote
I don't think your description is correct, at the moment you are solving a regular ODE. The solver is not aware that the last equation doesn't contain the derivative of x(5). And this needs to be specified with odeset 'Mass' option and other parameters
Yes, you need something like
M = diag([1 1 1 1 0]);
vopt = odeset ('InitialStep', 1e-3, 'AbsTol', 1e-8,'Mass',M);
I cannot test it since I do not have a working odepkg due to Octave 4.2 incompatibilities.

Lazy_Tom wrote
Out of curiosity, getting odest documentation (command 'doc odepkg' doesn't work on my windows release) and odepkg.pdf document has been a problem.
Actually "doc" has also never worked for me. Not sure, if this is a bug or an expected fail.

Lazy_Tom wrote
The pdf is mentioned in odeset help, odepkg help and website (http://octave.sourceforge.net/odepkg/overview.html) but to the best of my knowledge the actual location of it doesn't exist. In fact, after installing the odepkg I did a search odepkg.pdf on my octave directory with no results. What I've found though is an odepkg.texi file which is indeed the package manual but it's not a pdf.
The fact that odepkg.pdf is not correctly build, is a bug of the packaging script of odepkg. It seems the variable "TEXI2PDF = texi2pdf --clean" is not set in the Makefile. We will fix it in the next release.

Lazy_Tom wrote
This is how I finally found out about the 'Mass' option (yes there is a tiny example in odepkg_examples_dae but I don't think this is enough)
This is extensively explained in the pdf (e.g. page 14 in http://cosy.informatik.uni-bremen.de/sites/default/files/odepkg.pdf).
Reply | Threaded
Open this post in threaded view
|

Re: DAE daspk solver problem

Lazy_Tom
In reply to this post by Sebastian Schöps
Sebastian Schöps wrote
Lazy_Tom wrote
I don't know why DAE solvers index capabilities are not reported in Octave manual. Matlab at least tells you that you can go no further than index 1 (at least without doing some "magic"...)
You are right. Feel free to oben a bug report at savannah, maybe immediately suggesting a string that can be added to the help. However, I assume that dassl and daspk will be deprecated in some future version since we are moving to odeXY compatible solvers.
Hi Seb.

Oh, so there are other solvers that support full implicit ODE form? I see, any chance you can send me some links like you did for Radu 5?

About bug report, I don't know how to do it and I don't understand the "suggesting a string" (string = comment?), can you give me more details? Or, considering future deprecation

Sebastian Schöps wrote
Lazy_Tom wrote
about the paper you cited, yup, I'm aware of it and I agree, it is quite nice.
Well, at least I'll gave it a try, I think that for DAE problems I'll move to openModelica
I think that ode5r from odepkg should be able to solve index3 problems since its based on Radau 5 [1]. We hope to have a new release of odepkg quickly after Octave 4.2. I will keep in mind to include such information in the help.
I will probably give Radau a go, assuming I understood how to set up odeset :)
Cheers for the index info

Sebastian Schöps wrote
Regarding openModelica: yes, they do some clever tricks for higher index. As far as I recall it's based on dummy derivatives [2]. Such stuff is easier for modelica since it is in some sense a symbolic language. However, to my best knowledge, it is not really well suited for big problems. So, as usual, the best choice depends on what you want to do :)
Indeed. About index reduction, at the moment I can see it has three methods (except for the third, no idea what the others mean): Uode, Dynamic State Selection and Dummy Derivatives

BTW, nice link!

Cheers,

Marco
Reply | Threaded
Open this post in threaded view
|

Re: DAE daspk solver problem

Lazy_Tom
In reply to this post by Sebastian Schöps
Sebastian Schöps wrote
Lazy_Tom wrote
I don't think your description is correct, at the moment you are solving a regular ODE. The solver is not aware that the last equation doesn't contain the derivative of x(5). And this needs to be specified with odeset 'Mass' option and other parameters
Yes, you need something like
M = diag([1 1 1 1 0]);
vopt = odeset ('InitialStep', 1e-3, 'AbsTol', 1e-8,'Mass',M);
I cannot test it since I do not have a working odepkg due to Octave 4.2 incompatibilities.
I can give it a go. This is what I tried, it is based on the odepkg_examples_dae Demonstration 1 (which I know it works)

to make things simple m=1, g=10, pole length=1, both ode2r and ode5r has been tried.

function [vyd] = frobertson (vt, vy, varargin)
   vyd(1,1) =  vy(2);
   vyd(2,1) =  -vy(1) * vy(5);
   vyd(3,1) =  vy(4);
   vyd(4,1) =  -vy(3) * vy(5) - 10;
   vyd(5,1) =  vy(1)^2 + vy(3)^2 - 1;
endfunction


function [vmass] = fmass (vt, vy, varargin)
   vmass =  diag([1, 1, 1, 1, 0]);
endfunction

vopt = odeset ('Mass', @fmass, 'NormControl', 'on');
vsol = ode5r (@frobertson, [0, 17], [1, 0, 0, 0, 0], vopt);
plot (vsol.x, vsol.y);

the output was :

warning: Option "RelTol" not set, new value 1.0e-006 is used
warning: called from
    radau5_test_1_V04 at line 19 column 6
warning: Option "AbsTol" not set, new value 1.0e-006 is used
warning: Option "NormControl" will be ignored by this solver
warning: Option "InitialStep" not set, new value 1.0e-006 is used
warning: Option "MaxStep" not set, new value 1.4e+000 is used
warning: Option "Mass" only supports constant mass matrices M() and not M(t,y)
warning: Option "NewtonTol" not set, default value is used
warning: Option "MaxNewtonIterations" not set, default value 7 is used
 EXIT OF RADAU5 AT X=        0.2537E-02
  STEP SIZE T0O SMALL, H=   1.7937759160684990E-018
error: missing implementation
error: called from
    radau5_test_1_V04 at line 19 column 6

Unless I have made some typos, I suspect this might take more than expected to set up...


Sebastian Schöps wrote
Lazy_Tom wrote
This is how I finally found out about the 'Mass' option (yes there is a tiny example in odepkg_examples_dae but I don't think this is enough)
This is extensively explained in the pdf (e.g. page 14 in http://cosy.informatik.uni-bremen.de/sites/default/files/odepkg.pdf).
Cheers, I assume there are no changes in odeset from 0.8.2 to 0.8.5

Marco

Reply | Threaded
Open this post in threaded view
|

Re: DAE daspk solver problem

Sebastian Schöps
Lazy_Tom wrote
the output was :

warning: Option "RelTol" not set, new value 1.0e-006 is used
warning: called from
    radau5_test_1_V04 at line 19 column 6
warning: Option "AbsTol" not set, new value 1.0e-006 is used
warning: Option "NormControl" will be ignored by this solver
warning: Option "InitialStep" not set, new value 1.0e-006 is used
warning: Option "MaxStep" not set, new value 1.4e+000 is used
warning: Option "Mass" only supports constant mass matrices M() and not M(t,y)
Did you see the last line? The solver does not like function handles for the mass matrix. However, you can just use the matrix itself:

vmass =  diag([1, 1, 1, 1, 0]);
vopt = odeset ('Mass', vmass, 'NormControl', 'on');

I will try to look to your example in more detail next week.

Sebastian
Reply | Threaded
Open this post in threaded view
|

Re: DAE daspk solver problem

Lazy_Tom
Sebastian Schöps wrote
Lazy_Tom wrote
the output was :

warning: Option "RelTol" not set, new value 1.0e-006 is used
warning: called from
    radau5_test_1_V04 at line 19 column 6
warning: Option "AbsTol" not set, new value 1.0e-006 is used
warning: Option "NormControl" will be ignored by this solver
warning: Option "InitialStep" not set, new value 1.0e-006 is used
warning: Option "MaxStep" not set, new value 1.4e+000 is used
warning: Option "Mass" only supports constant mass matrices M() and not M(t,y)
Did you see the last line? The solver does not like function handles for the mass matrix. However, you can just use the matrix itself:

vmass =  diag([1, 1, 1, 1, 0]);
vopt = odeset ('Mass', vmass, 'NormControl', 'on');

I will try to look to your example in more detail next week.

Sebastian
The reasono why I didn't use the matrix and opted for the ready-made example is that the matrix didnt work.
In any case, I send you the code and the error, hope it helps

function [vyd] = frobertson (vt, vy, varargin)
   vyd(1,1) =  vy(2);
   vyd(2,1) =  -vy(1) * vy(5);
   vyd(3,1) =  vy(4);
   vyd(4,1) =  -vy(3) * vy(5) - 10;
   vyd(5,1) =  vy(1)^2 + vy(3)^2 - 1;
endfunction

vmass =  diag([1, 1, 1, 1, 0]);

vopt = odeset ('Mass', vmass, 'NormControl', 'on');
vsol = ode5r (@frobertson, [0, 17], [1, 0, 0, 0, 0], vopt);
plot (vsol.x, vsol.y);

output:

warning: Option "RelTol" not set, new value 1.0e-006 is used
warning: called from
    radau5_test_1_V05 at line 18 column 6
warning: Option "AbsTol" not set, new value 1.0e-006 is used
warning: Option "NormControl" will be ignored by this solver
warning: Option "InitialStep" not set, new value 1.0e-006 is used
warning: Option "MaxStep" not set, new value 1.4e+000 is used
warning: Option "NewtonTol" not set, default value is used
warning: Option "MaxNewtonIterations" not set, default value 7 is used
 EXIT OF RADAU5 AT X=        0.2428E-02
  STEP SIZE T0O SMALL, H=   2.1389523049476520E-018
error: missing implementation
error: called from radau5_test_1_V05 at line 18 column 6


Marco