

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


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.huberlin.de/~steffen/pub/introduction_to_daes_497.pdf (page 2 and 13).
Seb.


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


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


I was able to solve this using ode5r from odepkg. My test code is below:
function radau5PendulumDAE
vopt = odeset ('InitialStep', 1e3, 'AbsTol', 1e8); 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
_______________________________________________
Helpoctave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/helpoctave


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?


I did not try this simple test in daspk but I did try it in ida (with Capi). 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 index3 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'.
_______________________________________________
Helpoctave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/helpoctave


DASPK/DASKR are not expected to solve index 2+ problems like the index3 pendulum problem.
My comment is about problem structure (mass matrix vs fully implicit), sparsity and scaling.


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 matrixfree iterative solver
in DASPK from Octave.
c.
_______________________________________________
Helpoctave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/helpoctave


Bill Greene3 wrote
I was able to solve this using ode5r from odepkg. My test code is below:
function radau5PendulumDAE
vopt = odeset ('InitialStep', 1e3, 'AbsTol', 1e8);
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


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', 1e3, 'AbsTol', 1e8,'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.unibremen.de/sites/default/files/odepkg.pdf).


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


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', 1e3, 'AbsTol', 1e8,'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.0e006 is used
warning: called from
radau5_test_1_V04 at line 19 column 6
warning: Option "AbsTol" not set, new value 1.0e006 is used
warning: Option "NormControl" will be ignored by this solver
warning: Option "InitialStep" not set, new value 1.0e006 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.2537E02
STEP SIZE T0O SMALL, H= 1.7937759160684990E018
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...
Cheers, I assume there are no changes in odeset from 0.8.2 to 0.8.5
Marco


Lazy_Tom wrote
the output was :
warning: Option "RelTol" not set, new value 1.0e006 is used
warning: called from
radau5_test_1_V04 at line 19 column 6
warning: Option "AbsTol" not set, new value 1.0e006 is used
warning: Option "NormControl" will be ignored by this solver
warning: Option "InitialStep" not set, new value 1.0e006 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


Sebastian Schöps wrote
Lazy_Tom wrote
the output was :
warning: Option "RelTol" not set, new value 1.0e006 is used
warning: called from
radau5_test_1_V04 at line 19 column 6
warning: Option "AbsTol" not set, new value 1.0e006 is used
warning: Option "NormControl" will be ignored by this solver
warning: Option "InitialStep" not set, new value 1.0e006 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 readymade 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.0e006 is used
warning: called from
radau5_test_1_V05 at line 18 column 6
warning: Option "AbsTol" not set, new value 1.0e006 is used
warning: Option "NormControl" will be ignored by this solver
warning: Option "InitialStep" not set, new value 1.0e006 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.2428E02
STEP SIZE T0O SMALL, H= 2.1389523049476520E018
error: missing implementation
error: called from radau5_test_1_V05 at line 18 column 6
Marco

