# DAE daspk solver problem

15 messages
Open this post in threaded view
|

## DAE daspk solver problem

 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
Open this post in threaded view
|

## Re: DAE daspk solver problem

 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.
Open this post in threaded view
|

## Re: DAE daspk solver problem

 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
Open this post in threaded view
|

## Re: DAE daspk solver problem

 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
Open this post in threaded view
|

## Re: DAE daspk solver problem

 I was able to solve this using ode5r from odepkg. My test code is below:function radau5PendulumDAEvopt = 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));endfunction 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];endOn Fri, Oct 28, 2016 at 2:18 AM, 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. 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
Open this post in threaded view
|

## Re: DAE daspk solver problem

 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?
Open this post in threaded view
|

## Re: DAE daspk solver problem

 I did not try this simple test in daspk but I did try it in ida (with C-api).Even after carefullyinsuring that I had consistent initial conditions, ida could not advance to thefirst time step. So I am very doubtful that daspk will solve this example in itsindex-3 form.Frankly, I was quite surprised that radau5 was able to do so particularly sinceit doesn't appear to have an option to specify an initial y'. On Fri, Oct 28, 2016 at 1:17 PM, SunilShah 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
Open this post in threaded view
|

## Re: DAE daspk solver problem

 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.
Open this post in threaded view
|

## Re: DAE daspk solver problem

 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
Open this post in threaded view
|

## Re: DAE daspk solver problem

 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
Open this post in threaded view
|

## Re: DAE daspk solver problem

 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).
Open this post in threaded view
|

## Re: DAE daspk solver problem

 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
Open this post in threaded view
|

## Re: DAE daspk solver problem

 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