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

2 次查看(过去 30 天)
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 个评论
Christian Gärtner
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!

请先登录,再进行评论。

采纳的回答

darova
darova 2020-2-18
编辑:darova 2020-2-18
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 个评论
Christian Gärtner
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
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

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by