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

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!

[darova]

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]

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]

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

[Christian Gärtner]

thank you! that really helped me a lot.:)

best regards!

[darova]

you are welcome

[Christian Gärtner]

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

[darova]

• 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));

As i inderstood you correctly TL is kind of force (Newtons unit)

TL0-y(2)            % y(2) is force too (Newtons?)

So may i know what if Fext(units?)?

Also, you don't for loop here

[Christian Gärtner]

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]

• 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