Unable to perform assignment because the left and right sides have a different number of elements
2 次查看(过去 30 天)
显示 更早的评论
% 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 个评论
KSSV
2022-3-10
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.
回答(1 个)
Torsten
2022-3-10
编辑:Torsten
2022-3-15
% 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 个评论
Torsten
2022-3-15
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
2022-3-15
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)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Combustion and Turbomachinery 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!