Bendesarts

345 Reputation

10 Badges

10 years, 122 days

MaplePrimes Activity


These are replies submitted by Bendesarts

@tomleslie 

Thanks a lot for your very detailled and relevant help.

Concerning the points display, it comes from the output=listprocedure in the dsolve function. I attached this file in order to show you what it was done before to plot points on the ellipse GaitPlot_points.mw

May you see if you can correct the file GaitPlot5 file so as to enable the display of some points (at t = 0, 0.5, 1, 1.5 as in the file attached) ?

Thank you for your help

@ Dr. Venkat Subramanian

Thank you for you answer.

I agree that MapleSim have DAE solver.

But, are you sure that the NL oscillator is a DAE ?

For my me, it is a set of 8 coupled ODES but not a DAE. Indeed, there is not algebraic equations.

But, you probably have far more experiences than me on DAE.

@tomleslie 

OK. Thank you for you feedback. I try to follow your tips. I try to apply you method on the NL oscillator that I use in my robot model.

But, I think that slight modifications on the calculation of the vertical space vSpace should be conducted.

In order to facilitate the tuning of the display of the plots, I attach here a worksheet where I added comments so as to facilitate the reading.

GaitPlotTrajectory.mw

I hope you will find this project more interesting and I thank you again for your help.

 

 

@tomleslie 

Hello,

Thank you again for your help and relevant advices.

I'm very close to the result expected.

1) For the adaptation to the parameters modification

Thanks to the parameter vSpace, I almost obtain a "nice" result from the beginning. You can see that I need to add a coefficient with a value of 5 to get a good plot. I will be very interested if you see how I can adapt my plot without adding this value of 5 but with a more appropriate parameter. In order to help you to have a better understanding of my NL oscillator, I put the code at the end of this message. It will enable you not to make the magicien to understand my parameters. Concerning vSpace, the NL oscillator is design to make an ellipse with a heigh 2H and a length L. So, I think that vSpace:=2H should fit well. What do you think about that ?

2) Adding vertical lines 

Another interesting display could be to add vertical lines at the beginning of each red phase. In fact, the red phase correspond to the phase where the feet are in contact with the floor.

Here a example of the vertical lines that I would like to add :

Do you have ideas so as to add these vertical lines ?

3) Points on the ellipse

There is a element which does no longer work. It is the addition of points on the ellipse. I think it comes from the fact that I use now the dsolve with output=listprocedure. Do you have ideas so as to obtain the coordinates p0, p1, p2 and p3 in this case where res is obtained with this dsolve with output=listprocedure ?

K:=Matrix([<0, -1, -1, 1>,<-1, 0, 1, -1>,<1, -1, 0,-1>,<-1, 1, -1,0>]);
omega[sw]:=beta/(1-beta)*omega[st];
for i to 4
do
r[i]:=sqrt((u[i](t)/(L/2))^2+(v[i](t)/H)^2):
omega[i]:=omega[st]/(1+exp(b*v[i](t)))+omega[sw]/(1+exp(-b*v[i](t))):
Equ[i]:=diff(u[i](t),t)=Au*(1-r[i]^2)*u[i](t)+omega[i]*(L/2)/H*v[i](t):
Eqv[i]:=diff(v[i](t),t)=Av*(1-r[i]^2)*v[i](t)-omega[i]*H/(L/2)*u[i](t)+MatrixVectorMultiply(K,<seq(v[i](t),i=1..4)>)[i]:
EqSys[i]:=[Equ[i],Eqv[i]]:
end do:
EqSys[1];
paramsGeo:=L=0.080,H=0.020,beta=0.5,Vf=0.3;
omegaS:=eval(Pi*Vf/L, [paramsGeo]);
paramsCycle:=omega[st]=omegaS,Au=5,Av=50,b=100;
params:=paramsGeo,paramsCycle;
sys:=map(op,eval([seq(EqSys[i],i=1..4)],[params]));
ic:=[u[1](0)=-0.02, v[1](0)=0,u[2](0)=-0.015, v[2](0)=0,u[3](0)=0, v[3](0)=0.01,u[4](0)=0, v[4](0)=0.01];
Resolution
res:=dsolve([sys[],ic[]],numeric, maxfun = 0);
Initial boundaries calculations
tcalc:=8;
ic:=[seq(u[i](0)=eval(u[i](t), res(tcalc)),i=1..4),seq(v[i](0)=eval(v[i](t), res(tcalc)),i=1..4)];
Resolution2
res:=dsolve([sys[],ic[]],numeric, maxfun = 0,output=listprocedure);
Tracé
tmax:=10;

numpts:=100*tmax;
plots:-odeplot(res,[u[1](t),v[1](t)],t=0..tmax,scaling=constrained,numpoints=numpts);
plt1 := %:
plots:-odeplot(res,[u[2](t),v[2](t)],0..tmax,scaling=constrained,numpoints=numpts);
plots:-odeplot(res,[u[3](t),v[3](t)],0..tmax,scaling=constrained,numpoints=numpts);
plots:-odeplot(res,[u[4](t),v[4](t)],0..tmax,scaling=constrained,numpoints=numpts);

Points on ellipse
p0:=eval([u[1](t),v[1](t)], res(0))
p1:=eval([u[1](t),v[1](t)], res(0.5));
p2:=eval([u[1](t),v[1](t)], res(1));
p3:=eval([u[1](t),v[1](t)], res(1.5));
pointplot([p0, p1, p2,p3],color=[red,green,blue,pink], symbol=solidcircle, symbolsize=20); plt2 := %:
display([plt1, plt2]);
plots:-odeplot(res,[seq([t,v[i](t)+i/10],i=1..4)],0..tmax,thickness=2,numpoints=numpts);
vSpace:=eval(H, [paramsGeo]);
tmax:= 10:
numpts:=100*tmax:
TraceCourbe:= plots:-odeplot(res,[seq([t,v[i](t)+(i-1)*5*vSpace],i=1..4)],0..tmax,thickness=2, view=[0..5, -1*5*vSpace..4*5*vSpace], numpoints = numpts):
tmax:=5:
for j from 1 by 1 to 4 do
v:= rhs(res[5+j]);
la||j:= plot( [seq(`if`(v(t)>0,[t, vSpace*5*(j-1)],NULL),t=0..tmax,0.01)],
style=point,
color=blue,
symbol=solidcircle,
symbolsize=2
):
lb||j:= plot( [seq(`if`(v(t)<0,[t, vSpace*5*(j-1)],NULL), t=0..tmax,0.01)],
style=point,
color=red,
symbol=solidcircle,
symbolsize=2
):
od:

plots[display]([TraceCourbe, seq([la||j, lb||j][], j=1..4)],tickmarks=[default, [0=0, 0.1=0, 0.2=0, 0.3=0]],axes=boxed);

Thanks a lot for your help

@tomleslie 

Thank a lot for your help.

For 1) about the x-axis, it is now perfect with the option axes=boxed 

For 2) Your technique for finding the minimum and maximum are great. Thank you. However, I would have liked to display this max and min value in the graph in fact. 

3) I meet another problem. When I change the parameters of the NL oscillator, I need to adapt the term in bold in the following equation:

p1:= plots:-odeplot(res,[seq([t,v[i](t)+(i-1)*5],i=1..4)],0..tmax,thickness=2, view=[0..5, -1..16], numpoints = numpts):

Do you have ideas so as to have a more scalable (or flexible) way in case of modifications the parameters of the oscillator ? The idea is to adapt the plotting in case of curves with different amplitudes.

Thank you in advance for your ideas.

 

@tomleslie 

Your solution "Edit->Remove Output->From Worksheet" works perfect !

Thank you

@tomleslie 

Thanks a lot for your help.

It is almost perfect for me.

The only points remaining are :
- for the the blue/ red lines associated to v[1](t), the blue and red colors are not easy to see because the x-axis also on this line.

- Concerning the y-axis, the display is not so bad, but it is a pity that it is not possible to have a appropriate graduation by using this method. I could be interested to analyse for example the amplitude of a maximum which is a bit difficult here because the graduation of the y-axis is quiet limit.

 

@tomleslie 

Thanks a lot for you help.

I find the result not so bad and not too far from what I need.

Here the result obtained with you method :

1) About the code

First, I try to understand your code :

tmax:=5:
for j from 1 by 1 to 4 do
v:= rhs(res[5+j]);
la||j:= plot([seq(`if`(v(t)>0,[t, 5*j],NULL),t=0..tmax,0.01)], color=blue, style=point):
lb||j:= plot([seq(`if`(v(t)<0,[t, 5*j],NULL),t=0..tmax,0.01)], color=red, style=point):
od:
plots[display]([p1, seq([la||j, lb||j][], j=1..4)]);

May you explain me what the code res[...] means ? In fact, I usually use res() for asking a value but not with brackets "[]" (Nota : res:=dsolve(...))

2) About the gait diagrams

This result almost fits with my expectation.

a) Do you have ideas so as to make the blue and the red line smaller ?

b) Is it possible to plot this kind of gait diagrams by starting for the vertical axes each time by 0 like in this way :

 

Thanks a lot for you help

 This subject may interest you. That's why I put you informed.

@Carl Love 

@Preben Alsholm 

@Preben Alsholm 

@tomleslie 

@Carl Love 

Thank you for your support.

Concerning the remarks of @tomleslie : 

1) I think that you are right that behind MapleSim almost all the functions used are based or directly are Maple functions.

2) OK, may be. But the movement of the foots of the robots (elliptic trajectories) have directly a impact on the general trajectory of the robot's body. And after a moment of simulation (around 25 s), i see the robot which don't follow a straight line and I guess that it comes from the elliptic trajectories of the foots which are no longer purely elliptic.

Concerning the remarks of @Preben Alsholm :

1) I have indeed tried to use Maple first and besides the code used in MapleSim comes from my Maple's code. The only difference is the fact, in MapleSim, I no longer use the dsolve and odeplot functions but this code :

sModelica := MapleSim:-Tools:-MapleToModelica(
[seq(EqSys[i],i=1..4)],
'initial'=ic,
'class_name' = "ClassTrajectoire",
'comment' = "Trajectoire utilisée pour le mouvement des pieds",
'parameters' = params,
'connections' = cons,
'unfilled_pins' = {"u1_out", "v1_out","u2_out","v2_out","u3_out","v3_out","u4_out","v4_out"}
):

and the solver of MapleSim.

2) In Maple, I obtain limit cycles I believe. So, normally, I should obtain also limit cycles in MapleSim

Don't hesitate to let me know if you have new ideas for troubleshooting my issue of "drift" for the solution of this NL oscillator but in MapleSim.

Thanks a lot for your help 

@Preben Alsholm 

@tomleslie 

@Carl Love 

Thank for your ideas on the solution of the non linear oscillator for obtaining the elliptic trajectory.

I'm still confused about my results obtained with MapleSim.

In fact, I totally agree with you about the issue of visualisation with the plot function when I don't ask enough points with Maple.

However, in MapleSim, I follow the tip of @Preben Alsholm by using the initial conditions obtained after 8s. To doing this, I launch one time the oscillator. I retrieve the values of u[i], v[i] for t=8s and I use these values as initial conditions for the next run. So, I use directly good initial conditions for solving the differential system.

Here the result I obtain with MapleSim:

 

It thus seems that I observe a "drift" on the movement of the robot's model. Indeed, my robot has a perfect trajectory until around 25s of simulation. Besides, until 25s, the trajectory obtained with the NL oscillator is purely elliptic. But, after this time, it begins to diverge slowly.

It may come from the solver used by MapleSim.

Do you have some ideas so as to tune the oscillator or the solver used in MapleSim so as to keep the periodic trajectory for a long time of simulation ?

I can send you the different solvers which can be used with MapleSim.

Thanks a lot for your help.

@tomleslie 

@Carl Love 

@Preben Alsholm 

Thanks a lot for your help.

You have given me a method if I want to obtain a nice visualization of the solution of my oscillator in Maple. This solution is based on a appropriated choice of numpoints with the plot function.

However, I'm wondering if there is only a "visual" mistake and not also some drifts during the simulation.

The reason for this guess is the fact that I also obtain some deviations when I generate my elliptic trajectory with MapleSim. My MapleSim code is based on my maple code. I use directly the equations from my maple code into MapleSim. However, the solver is different in the sense that I didn't use the dsolve function from maple but the solver of mapleSim.

My MapleSim code is :

sModelica := MapleSim:-Tools:-MapleToModelica(
[seq(EqSys[i],i=1..4)],
'initial'=ic,
'class_name' = "ClassTrajectoire",
'comment' = "Trajectoire utilisée pour le mouvement des pieds",
'parameters' = params,
'connections' = cons,
'unfilled_pins' = {"u1_out", "v1_out","u2_out","v2_out","u3_out","v3_out","u4_out","v4_out"}
):

It enables to create a dynamical system which can be solved by MapleSim solver

The settings of the MapleSim solver is :

The result obtained with MapleSim is :

 

1) In maple :

If there is really a drift during the resolution made by the dsolve function, do you have ideas to reduce this drift so that I can, in my case, obtain a nice ellipse for my generated trajectory?

In others words, in the previous posts, we have tried to better tune the plot function. But, before trying to correct the visualization, is it possible to correct the deviations at the "dsolve" function level in maple ?

2) In mapleSim:

Do you have ideas so as to tune properly the solver so as to avoid this drift from the elliptic trajectory ?

Thanks a lot for your help.

@tomleslie 

@Preben Alsholm 

Thanks a lot for your feedback.

However, for my application (robotics), despite the fact that the solution of the oscillator can be periodic (in some cases), I do believe that I need to obtain the solution for the duration of my simulation. I will argue with 3 reasons:

1) The oscillators which I want to use contain 2 half period (one for the stance phase and one for the support phase). We can see the definitions of 2 half periods in the second example that I gave.

2) I have the simulations conducted also with matlab. In matlab, it seems that there is no problem to ask a big tmax. Why cannot I do the same with Maple? I hope that I can tune dsolve in Maple in order to have a solution from a larger range of t.

3) These kind of oscillators are very useful in robotics because it is possible to add a sensory feedback signal (ui) in the oscillator like in the following equations :

Consequently, the oscillator is slightly pertubed to take into account some information from environment but comes back after a moment to its limit cycle. This is the kind of oscillators that I'm trying to tune in Maple.

So, I hope that my elements bring you more useful information and permit a mutual understanding. 

If yes, i will be very interested by you experience and feedback so as to tune "dsolve" (or another function) for having a simulation on a large range of time.

Do you have other ideas so as to obtain a solution for a large range of time and keeping a nice elliptical plot ?

Thank you for your help.

@Preben Alsholm 

Thank you for your help.

But, in fact, I use this non linear coupled oscillator for a robotic application. In my application, this oscillator enables me to calculate the trajectory of the foot of my robot.

Consequently, I need to make my simulation from t=0 to t=300.

So, my questions :

1) Do you have ideas so as to tune numpoints correctly ? I would like to build a rule of thumb so as to test directly a more appropriated number of point (numpoints).

2) As it is for the generation of trajectory for a robot, I'm quiet interesting to start the plotting when I have a good phase shift between my plots. Indeed, i put bad initial conditions but after 2 or 3s, the oscillator reachs its limit cycle.

May I do with this syntax : res:=dsolve({EqSys[],ic[]},numeric,range=tmin..tmax); with the definition of a tmin value ?

@tomleslie 

Perfect, thanks a lot for your very constructive and accurate help.

With this non linear oscillator, tuning the numpoints(with numpoints=10*tmax) enables to keep the solution with the limit cycle.

However, I tried to used the same tuning (that is to say to use the same number of points numpoints) in a more complex coupled oscillator and this time it doesn't work.

My code is the following :

K:=Matrix([<0, -1, 1, -1>,<-1, 0, -1, 1>,<-1, 1, 0,-1>,<1, -1, -1,0>]);

for i to 4
do
r[i]:=sqrt((u[i](t))^2+(v[i](t))^2):
omega[i]:=omega[sw]/(1+exp(b*v[i](t)))+omega[st]/(1+exp(-b*v[i](t))):
Equ[i]:=diff(u[i](t),t)=Au*(1-r[i]^2)*u[i](t)-omega[i]*v[i](t):
Eqv[i]:=diff(v[i](t),t)=Av*(1-r[i]^2)*v[i](t)+omega[i]*u[i](t)+MatrixVectorMultiply(K,<seq(v[i](t),i=1..4)>)[i]:
EqSys[i]:=[Equ[i],Eqv[i]]:
end do:

paramsCycle:=omega[st]=2*Pi,omega[sw]=4*2*Pi,Au=5,Av=5,b=100;
params:=paramsCycle;

sys:=map(op,eval([seq(EqSys[i],i=1..4)],[params]));
tmax:= 100:
ic:=[u[1](0)=0, v[1](0)=0,u[2](0)=0, v[2](0)=-0.1,u[3](0)=0, v[3](0)=0.1,u[4](0)=0, v[4](0)=0.1];

res:=dsolve([sys[],ic[]],numeric, maxfun = 0, output = listprocedure):

plots:-odeplot(res,[u[1](t),v[1](t)],0..tmax, scaling = constrained, numpoints = 10*tmax);
plots:-odeplot(res,[seq([t,v[i](t)+i/10],i=1..4)],0..tmax,thickness=2);

However, it seems to work with numpoints = 100*tmax

Do you have ideas so as to tune numpoints correctly ? I would like to build a rule of thumb so as to test directly a more appropriated number of point (numpoints).

Thank you in advance for your feedback.

 

On the result of the non linear oscillator, I have his red arrows. They should be correspond to the general solution. Which option can I set in order to keep only my limit cycle?

Thank your for your help.

First 6 7 8 9 10 11 12 Last Page 8 of 16