Code only works as expected with 1 value, else, "FunctionLine update: Unable to convert expression into double array" error
์กฐํ ์: 2 (์ต๊ทผ 30์ผ)
์ด์ ๋๊ธ ํ์
I need specific data points as output to feed into another project. I got to the point where i get the correct output when B (correction factor) is set to 0.8. The correct output should be 0 until tau in seconds, than 2 triangle shapes, as shown:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1609671/image.png)
The relevant code is as follows:
%Initializing the values for computation
R = 0.05;
H = 4.0;
K = 0.95;
Fh = 0.3;
Tr = 8.0;
D = 1.0;
Ps = 0.3; %Initializing the P-step value
%Initializing t and s
syms t s
%Calculating ๐๐
wn2 = (D*R+K)/(2*H*R*Tr);
wn = sqrt(wn2);
%Calculating ๐
c = (((2*H*R)+(((D*R)+(K*Fh))*Tr))/(2*(D*R+K)))*wn ;
%Calculating alpha
a = sqrt((1-2*Tr*c*wn+(Tr^2)*(wn^2))/(1-(c^2)));
%Calculating ๐r
wr = wn*(sqrt(1-(c^2)));
%Calculating the angle (In Radians)
phi = atan((wr*Tr)/(1-c*wn*Tr))-atan(sqrt(1-(c^2))/(-c));
%frequency response
ft = vpa(1 - ((R*Ps)/(D*R+K))*(1+a*exp(-c*wn*t)*sin(wr*t+phi))); %THIS SEEMS TO BE THE BIG OFFENDER
%Defining the frequency Nadir
fs = 1 - ((R*Ps)/(D*R+K));
tn = (1/wr)*atan((wr*Tr)/(c*wn*Tr-1));
fn = 1 - ((R*Ps)/(D*R+K))*(1+a*exp(-c*wn*tn)*sin(wr*tn+phi));
%Declare the correction Factor
B = 0.8;
Dtr = fs - fn;
%The new frequency Nadir with the correction Factor
fnew = B*Dtr + fn;
%Declare and compute Tau and Tau2
sympref('HeavisideAtOrigin',1);
tau = vpasolve(ft-(B*Dtr+fn)==0,t)
tau2 = vpasolve(ft-(B*Dtr+fn)==0,t,tau+tn)
%Define the function Delta Frequency
dft = (fnew - ft)*(heaviside(t-tau) - heaviside(t-tau2));
%final equation to get power injection
dFs = laplace(dft)
GSFR = ((R*(wn^2))/(D*R+K))*((1+Tr*s)/((s^2)+2*c*wn*s+(wn^2)));
dPs = dFs/GSFR;
dpt = ilaplace(dPs)
fplot(dpt, [0 20])
xlim([0 20])
ylim([-0.05 0.3])
I used vpa on the calculation for dft because without it, even 0.8 does not work.
In theory, I should be able to adjust the value for B between 0 and 1, and the value for Ps between 0 and 1. But if i don't use B=0.8 and Ps=0.3, i get the following error on the final fplot call:
Warning: Error updating FunctionLine.
The following error was reported evaluating the function in FunctionLine update: Unable to convert expression into double array.
I don't know what I did wrong. is it too much data for matlab? complexity of the equation too high? (I am modelling off a published paper, so the equations are not mine, the implementation however is)
Edit: I found some other pairings that work. B=0.5 and Ps = 0.2 or 0.15 works, and B=0.2 and Ps = 0.1 works.
Edit 2: So laplace() is unable to resolve t, as shown, and it doesnt work in Matlab 2023, with the numbers where it does work in Matlab 2019b.
Edit 3: Still actively debugging. it seems that the vpa is required for the function ft. It is also apparent that after the laplace() there remains instances of the variable t. I believe that narrows the window of problems to that function, ft, and why it is not properly converting to the laplace domain.
๋๊ธ ์: 4
Star Strider
2024๋
2์ 6์ผ
Looking at the expression for โdptโ, the e argument (exponent) is still a function of s, and s remains in the expressions of
and
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1609701/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1609706/image.png)
See how the authors of the paper managed to invert it.
๋ต๋ณ (1๊ฐ)
Pramil
2024๋
2์ 16์ผ
Hey, as you mentioned, the problem is indeed happening because of the โlaplaceโ function which is unable to convert the โdftโ entirely, due to which you are getting the error.
Here is a link to the MATLAB answer which resolves a similar issue :
You can define โtโ as a real number which is true in context of Laplace transformation and use the โsimplify()โ function to get the desired result.
Here is the MATLAB code that works in both, R2023b and R2019b:
%Initializing the values for computation
R = 0.05;
H = 4.0;
K = 0.95;
Fh = 0.3;
Tr = 8.0;
D = 1.0;
Ps = 0.3; %Initializing the P-step value
%Initializing t and s
syms t real;
syms s;
%Calculating ๐๐
wn2 = (D*R+K)/(2*H*R*Tr);
wn = sqrt(wn2);
%Calculating ๐
c = (((2*H*R)+(((D*R)+(K*Fh))*Tr))/(2*(D*R+K)))*wn ;
%Calculating alpha
a = sqrt((1-2*Tr*c*wn+(Tr^2)*(wn^2))/(1-(c^2)));
%Calculating ๐r
wr = wn*(sqrt(1-(c^2)));
%Calculating the angle (In Radians)
phi = atan((wr*Tr)/(1-c*wn*Tr))-atan(sqrt(1-(c^2))/(-c));
%frequency response
ft = vpa(1 - ((R*Ps)/(D*R+K))*(1+a*exp(-c*wn*t)*sin(wr*t+phi))); %THIS SEEMS TO BE THE BIG OFFENDER
%Defining the frequency Nadir
fs = 1 - ((R*Ps)/(D*R+K));
tn = (1/wr)*atan((wr*Tr)/(c*wn*Tr-1));
fn = 1 - ((R*Ps)/(D*R+K))*(1+a*exp(-c*wn*tn)*sin(wr*tn+phi));
%Declare the correction Factor
B = 0.4;
Dtr = fs - fn;
%The new frequency Nadir with the correction Factor
fnew = B*Dtr + fn;
%Declare and compute Tau and Tau2
sympref('HeavisideAtOrigin',1);
tau = vpasolve(ft-(B*Dtr+fn)==0,t);
tau2 = vpasolve(ft-(B*Dtr+fn)==0,t,tau+tn);
%Define the function Delta Frequency
dft = (fnew - ft)*(heaviside(t-tau) - heaviside(t-tau2));
%final equation to get power injection
dFs = laplace(dft);
dFs = simplify(dFs);
GSFR = ((R*(wn^2))/(D*R+K))*((1+Tr*s)/((s^2)+2*c*wn*s+(wn^2)));
dPs = dFs/GSFR;
dpt = ilaplace(dPs);
fplot(dpt, [0 20])
xlim([0 20])
ylim([-0.05 0.3])
๋๊ธ ์: 0
์ฐธ๊ณ ํญ๋ชฉ
์นดํ ๊ณ ๋ฆฌ
Help Center ๋ฐ File Exchange์์ Calculus์ ๋ํด ์์ธํ ์์๋ณด๊ธฐ
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!