Code only works as expected with 1 value, else, "FunctionLine update: Unable to convert expression into double array" error

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:

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.

[Dyuman Joshi]

MATLAB is unable to resolve the ilaplace() of the provided input to an explicit expression of t (see the edit above).

Subsequently, fplot() is unable to process the input provided to get a number that is to be used for plotting (which is what the error states as well)

You could try defining the values as symbolic numbers via sym and check if that works.

[Chaim Rui van de Ven]

did I not do that by declaring s and t above in

%Initializing t and s

syms t s

if not, how would i go about defining them? and is it just by chance that ilaplace can do it for some sets of values but not others?

Edit: Did Not Work. using sym made no difference to declaring with syms. Appreciate the help. I did notice also, the inability to resolve to an explicit expression of t first occurs with laplace(), not ilaplace()

[Star Strider]

Looking at the expression for ‘dpt’, the e argument (exponent) is still a function of s, and s remains in the expressions of and .

See how the authors of the paper managed to invert it.

[Chaim Rui van de Ven]

unfortunately, they don't expand on that. they merely state that Delta Frequency in the Laplace domain (dFs) divided by GSFR gives Delta Power in the laplace domain, and "By computing the inverse Laplace transform of ∆𝑃(𝑠), one can easily derive ∆𝑝(𝑡)."

im confused because, in theory, there should be no t variables in the equation for dFs, and there should be no s variables in the equation for dpt.

Answer 1

[Pramil]

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 :

https://www.mathworks.com/matlabcentral/answers/100244-why-does-the-laplace-function-not-work-for-certain-exponential-expressions-in-symbolic-math-toolbox

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])