Unable to perform assignment because the left and right sides have a different number of elements
1 view (last 30 days)
Show older comments
% Input Data Parameters
h=1; % Step size [s]
tf=50; % Solve Space of t=[0,tf] [s]
Ta=297;
P=[1197, 1.8, 15.8];
I = 2;
C=4;
% Simulation time interval
t=0:h:tf;
% Coefficient of irreversible and reversible terms
A=[-1298*10^(-3),4038*10^(-3),-4674*10^(-3),2500*10^(-3),-618.6*10^(-3),128.5*10^(-3)];
B=[8.814*10^(-3),23.01*10^(-3),-19.09*10^(-3),4.099*10^(-3),1.011*10^(-3),-0.1937*10^(-3)];
% Intial conditions
Ts=zeros(1,length(t));
Ts(1)=Ta;
SOC=zeros(1,length(t));
Relec = zeros(1,length(t));
dUOC = zeros(1,length(t));
Qgen = zeros(1,length(t));
SOC(1)=1;
% ODE equation to solve
f=@(t,Ts) (Ta-Ts)/(P(1).*(P(2)+P(3)))+(Qgen.*P(3))/(P(1).*(P(2)+P(3)));
% RK4 Loop
for i=1:(length(t)-1)
SOC(i+1) = SOC(i)-(1/C)*I*((t(i+1)-t(i))/3600);
Relec = A(1).*(SOC.^5)+A(2).*(SOC.^4)+A(3).*(SOC.^3)+A(4).*(SOC.^2)+A(5).*(SOC)+A(6);
dUOC = B(1).*(SOC.^5)+B(2).*(SOC.^4)+B(3).*(SOC.^3)+B(4).*(SOC.^2)+B(5).*(SOC)+B(6);
Qgen =(I^2) .* Relec-I.*dUOC;
k1=f(t(i),Ts(i));
k2=f(t(i)+0.5*h,Ts(i)+0.5*k1*h);
k3=f(t(i)+0.5*h,Ts(i)+0.5*k2*h);
k4=f(t(i)+h,Ts(i)+k3*h);
Ts(i+1)=Ts(i)+(h/6)*(k1+2*k2+2*k3+k4);
1 Comment
KSSV
on 10 Mar 2022
You need to rethink on your code. This function:
f=@(t,Ts) (Ta-Ts)/(P(1).*(P(2)+P(3)))+(Qgen.*P(3))/(P(1).*(P(2)+P(3)));
Gives you a array as output. It should give you a scalar as output. Check this function.
Answers (1)
Torsten
on 10 Mar 2022
Edited: Torsten
on 15 Mar 2022
% Input Data Parameters
h=1; % Step size [s]
tf=50; % Solve Space of t=[0,tf] [s]
Ta=297;
P=[1197, 1.8, 15.8];
I = 2;
C=4;
% Simulation time interval
t=0:h:tf;
% Coefficient of irreversible and reversible terms
A=[-1298*10^(-3),4038*10^(-3),-4674*10^(-3),2500*10^(-3),-618.6*10^(-3),128.5*10^(-3)];
B=[8.814*10^(-3),23.01*10^(-3),-19.09*10^(-3),4.099*10^(-3),1.011*10^(-3),-0.1937*10^(-3)];
% Intial conditions
Ts=zeros(1,length(t));
Ts(1)=Ta;
soc = zeros(1,length(t));
Relec = zeros(1,length(t)-1);
dUOC = zeros(1,length(t)-1);
Qgen = zeros(1,length(t)-1);
soc(1)=1;
% ODE equation to solve
f=@(t,T,Q) (Ta-T)/(P(1).*(P(2)+P(3)))+(Q.*P(3))/(P(1).*(P(2)+P(3)));
for i=1:(length(t)-1)
SOC = soc(i);
Relec(i) = A(1).*(SOC.^5)+A(2).*(SOC.^4)+A(3).*(SOC.^3)+A(4).*(SOC.^2)+A(5).*(SOC)+A(6);
dUOC(i) = B(1).*(SOC.^5)+B(2).*(SOC.^4)+B(3).*(SOC.^3)+B(4).*(SOC.^2)+B(5).*(SOC)+B(6);
Qgen(i) =(I^2) .* Relec(i)-I.*dUOC(i);
k1=f(t(i),Ts(i),Qgen(i));
k2=f(t(i)+0.5*h,Ts(i)+0.5*k1*h,Qgen(i));
k3=f(t(i)+0.5*h,Ts(i)+0.5*k2*h,Qgen(i));
k4=f(t(i)+h,Ts(i)+k3*h,Qgen(i));
Ts(i+1)=Ts(i)+(h/6)*(k1+2*k2+2*k3+k4);
soc(i+1) = soc(i)-(1/C)*I*((t(i+1)-t(i))/3600);
end
plot(t,Ts)
end
Usually, you would also have to modify Qgen for each call of the function f depending on time.
For simplicity, I assumed Qgen to be constant on t(i) <= t <= t(i+1)
3 Comments
Torsten
on 15 Mar 2022
Qgen(i) =(I^2) .* Relec(i)-I.*dUOC(i);
If Qgen is heat generation, you can set this term (which is heat generation at time t(i), I guess) to a value of your choice.
Jan
on 15 Mar 2022
As usual I mention, that -4674*10^(-3) is an expensive power operation, while -4674e-3 is a cheap constant. I've learned programming on a ZX81 with 1kB RAM and avoiding EXP is my destiny.
A(1).*(SOC.^5)+A(2).*(SOC.^4)+A(3).*(SOC.^3)+A(4).*(SOC.^2)+A(5).*(SOC)+A(6)
% Horner scheme:
SOC * (SOC * (SOC * (SOC * (A(1) * (SOC) + A(2)) + A(3)) + A(4)) + A(5)) + A(6)
See Also
Categories
Find more on Combustion and Turbomachinery in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!