How can I solve a system of differential equations including a bounded proportional controller?

Question

Christian Gärtner 2020-2-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/506148-how-can-i-solve-a-system-of-differential-equations-including-a-bounded-proportional-controller

评论： darova 2020-3-12

采纳的回答： darova

在 MATLAB Online 中打开

Hello there,

I am kindly asking for help because I got stuck on this issue.

What I want to do is solve a set of differential equations which rely on a controller. The set of equations is as follows:

where M and F are constants. C(t) is the bounded proportional controller on which the system relies on, and which is given by the following:

where TL is a 2007x1 double vector. LD and LR is a constant.

I have tried to solve it by using the ode 45 function and by modelling the controller with an if-else statement. The script I have been using is as follows:

time=daten(8:end,2);
time_min=time/60;
T0=daten(8,2);
T0_min=T0/60;
tEnd=daten(end,2);
tEnd_min=tEnd/60;
Fext=daten(8:end,3)*100;    %_deltoideus_anterior_part1 in %
Fext_statisch=daten(8,3)*100; 
MVC=100;
TL=Fext/MVC;                 %Target Load in % of MVC
timerange=[T0_min tEnd_min];
IC=[MVC;0;0;0];
[t,y]=ode45(@(t,y) freylaw(t,y,time_min,TL),timerange,IC);
figure
plot(t,y(:,1));
hold on
plot(t,y(:,2));
hold on
plot(t,y(:,3));
xlabel('time[min]');
ylabel('Force of various compartments [N]');
legend('MR','MA','MF')

and the function I am calling in the ode45 statement is as follows:

function y=freylaw(t,y,time_min,TL)
L=10;                       %model-specific factor (LD=LR=5)
r=15;                       %rest-recovery-multiplier according to Looft et al.(2018)
F=0.00912;                  %fatigue factor according to Frey Law (2012)
R=0.00094;                   %recovery factor according to Frey Law (2012)
y(1)=-y(4)+(R+r)*y(3);      %Mr the constant r is introduced here according to Looft et al.(2018) and is an intended variation of the original equation
y(2)=y(4)-F*y(2);           %Ma
y(3)=F*y(2)-R*y(3);         %Mf
    if (y(2)<TL)&(y(3)>(TL-y(2)))
    y(4)=L*(TL-y(2));       %if Ma<TL and Mr>(TL-Ma)
    elseif (y(2)<TL)&(y(1)<(TL-y(2)))
    y(4)=L*y(1);            %if Ma<TL and Mr<(TL-Ma)
    else
    y(4)=L*(TL-y(2));       %if Ma>=TL
    end
end

The problem I am encountering is, that my result is just my initial condition stretched over the whole time range. But actually Ma should be a varying curve, and Mf should increase, and Mr should decrease. I guess that there is some issue with the if-statement I am using to model the bounded controller. But I am unsure on what to change and I am hoping some of you might give me some great tips! Thanks a lot!

2 个评论
显示无隐藏无

darova 2020-2-18

Can you please explain more:

what is TL and how it changes with time?
why do you use y(4) as C(t)?

Christian Gärtner 2020-2-18

编辑：Christian Gärtner 2020-2-18

Of course!

TL refers to the target load, which actually is the muscle activity a specific muscle (deltoid muscle in this case) is subject to while perfoming several movements. It is a value ranging from 0-1 according to 0-100%. you can see how it changes in the following plot:

I was thinking I should use y(4) instead of C(t) because I don't have a Value for C(t) until I have solved the differential equations. And I can't solve the differential equations until I have something I can use for C(t). So I thought making C(t) part of the system of differential equations will lead to desired result, though I am very unsure if that is the right way to do it.

Thank you for taking the time to think about the matter!

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

darova 2020-2-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/506148-how-can-i-solve-a-system-of-differential-equations-including-a-bounded-proportional-controller#answer_416144

编辑：darova 2020-2-18

在 MATLAB Online 中打开

Interpolate TL to choose correct value

Don't use y(4) for Ct

function dy=freylaw(t,y,time_min,TL)
L=10;                       %model-specific factor (LD=LR=5)
r=15;                       %rest-recovery-multiplier according to Looft et al.(2018)
F=0.00912;                  %fatigue factor according to Frey Law (2012)
R=0.00094;                   %recovery factor according to Frey Law (2012)
tt = linspace(0,1.6*60,length(TL));     % time array in seconds
TL0 = interp1(tt,TL,t);                 % choose current TL
if (y(2)<TL0)&&(y(1)>(TL0-y(2)))
    Ct=L*(TL0-y(2));       %if Ma<TL and Mr>(TL-Ma)
elseif (y(2)<TL0)&&(y(1)<(TL0-y(2)))
    Ct=L*y(1);            %if Ma<TL and Mr<(TL-Ma)
else
    Ct=L*(TL0-y(2));       %if Ma>=TL
end
dy(1,1)=-Ct+(R+r)*y(3);      %Mr the constant r is introduced here according to Looft et al.(2018) and is an intended variation of the original equation
dy(2,1)=Ct-F*y(2);           %Ma
dy(3,1)=F*y(2)-R*y(3);         %Mf
end

Small mistake

Edited: dy name of result variable

6 个评论
显示 4更早的评论隐藏 4更早的评论

Christian Gärtner 2020-3-11

在 MATLAB Online 中打开

Hello darova,

I kindly ask for your help again as I am trying improve the system of differential equation, which you have been helping me with.

I want to use an other variable TL in one of the following two ways (doesn`t matter which one):

TL=Fext/(100-y(3));

or

TL=Fext/(y(1)+y(2));

I have tried to do it this way in my function statement:

function dy=freylaw(t,y,time_min,Fext)
L=10;                                   %Modellspezifischer Faktor
r=15;                                   %rest-recovery-multiplier nach 2018#Looft et al.
F=0.00912;                              %Ermüdungsfaktor nach Frey Law (2012)
R=0.00094;                              %Erholungsfaktor nach Frey Law (2012)
%MVC=100;
for i=0:Fext(end)
MVC=100-y(3);
TL=Fext/MVC;
%TL=((Fext)/(100-y(3)));
%TL=((Fext)/(y(2)+y(1)));
TL0 = interp1(time_min,TL,t);           % choose current TL
    
if (y(2)<TL0)&&(y(1)>(TL0-y(2)))
    Ct=L*(TL0-y(2));                    %if Ma<TL and Mr>(TL-Ma)
elseif (y(2)<TL0)&&(y(1)<(TL0-y(2)))
    Ct=L*y(1);                          %if Ma<TL and Mr<(TL-Ma)
else
    Ct=L*(TL0-y(2));                    %if Ma>=TL
end
dy(1,1)=-Ct+(R+r)*y(3);                 %Mr the constant r is introduced here according to Looft et al.(2018) and is an intended variation of the original equation
dy(2,1)=Ct-F*y(2);                      %Ma
dy(3,1)=F*y(2)-R*y(3);                  %Mf
end
end

But this version and also the two outcommented lines of TL=Fext/(100-y(3)) or TL=Fext/(y(2)+y(1)) dont lead to reasonable results. Do you have any idea, what I am doing wrong?

Thanks in advance and best regards,

Christian

Christian Gärtner 2020-3-12

TL and Fext both have Newtons as unit and are a kind of force.

To go into detail:

Fext is the external load which a muscle has to handle while contracting
TL= Fext/MVC with MVC=maximum Force a muscle can output=100-y(2)
y(2)=part of the muscle which is becoming fatigued over the course of time (this leads to a decrease in MVC over time --> MVC=100-y(2))

Does that make sense to you?

The literature I am using for this is the following if you would like to take a closer look:

Xia & Frey Law (2008): theoretical approach for modeling peripheral muscle fatigue and recovery

thanks for helping!

darova 2020-3-12

Does that make sense to you?

Not really actually

(dimensionless value)

You can operate only with the same dimensions

MVC=100-y(2)) that means that y(2) and MVC are the same dimension (Newtons)

You can't substract values/number of different dimensions

You code looks ok. The only comment i have: check dimensions

请先登录，再进行评论。

How can I solve a system of differential equations including a bounded proportional controller?

2 个评论
显示无隐藏无

采纳的回答

6 个评论
显示 4更早的评论隐藏 4更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

How can I solve a system of differential equations including a bounded proportional controller?

2 个评论 显示 无隐藏 无

采纳的回答

6 个评论 显示 4更早的评论隐藏 4更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

2 个评论
显示无隐藏无

6 个评论
显示 4更早的评论隐藏 4更早的评论