Unable to perform assignment because the left and right sides have a different number of elements.

3 views (last 30 days)
% differential equation (SWING EQUATION)
%d^2(y)/d^2(t) = (-a*D)/(2*H)*dy/dt + (a*P)/(2*H) -(G*B*a)/(2*h*z)*sin(y-c);
%%MATLAB CODE%%
clc;
clear all;
N=100000;
T=linspace(0,100,N);
H=T(2)-T(1);
c=0.1; % theta B
H1=2.37; %Inertia
D= 0.002: 0.016; %damping term
P=0.85; %power
a=1; %constant angular velocity
G=1; %voltage of bus
f=60; %excitation frequesncy
s1=1; %psi v
B=1 + 0.1*cos(f*T + s1); %voltage of bus
z=0.645; %transient reactance
F=@(X2)X2;
G=@(X1,X2,T)-(a.*D)./(2.*H1).*X2 + (a.*P)./(2.*H1) -(G.*B.*a)./(2.*H1.*z).*((X1 - (((X1).^3)./6)).*cos(c) -(1-(((X1).^2)./2).*sin(c)));
X1=zeros(1,N);
X2=zeros(1,N);
X1(1)=2;
X2(1)=0;
for i=1:N-1
K1=F(X2(i));
L1=G(X1(i),X2(i),T(i));
K2=F(X2(i)+((H/2)*L1));
L2=G(X1(i)+((H/2)*K1),X2(i)+((H/2)*L1),T(i));
K3=F(X2(i)+((H/2)*L2));
L3=G(X1(i)+((H/2)*K2),X2(i)+((H/2)*L2),T(i));
K4=F(X2(i)+((H)*L3));
L4=G(X1(i)+((H)*K3),X2(i)+((H)*L3),T(i));
X1(i+1)=X1(i)+(H/6)*(K1+2*K2+2*K3+K4);
X2(i+1)=X2(i)+(H/6)*(L1+2*L2+2*L3+L4);
end
figure(1)
plot(T,X1)
hold on
plot(T,X2,'r');
legend('x_1(t)','x_2(t)');
xlabel('Normalized Time');ylabel('States of the system')
figure(2)
plot(X1,X2)
xlabel('x_1');ylabel('x_2')
Results:
Unable to perform assignment because the left and right sides have a different number of elements.
Error in TRIAL (line 31)
X1(i+1)=X1(i)+(H/6)*(K1+2*K2+2*K3+K4);
Can someone please help me with this code? I do not know where I went wrong. It is the swing eqaution and i need to sove using rungekutta method. The output i get is also methioned above!
  7 Comments
VBBV
VBBV on 23 Jun 2022
Edited: VBBV on 23 Jun 2022
B=1 + 0.1*cos(f*T + s1); %voltage of bus
In your definition it is vector with size of 1x100000. Try with a scale as below , it should work.
t = 1;% e.g.
B=1 + 0.1*cos(f*t + s1); %voltage of bus
KSSV
KSSV on 23 Jun 2022
I suspect B should be having a time difference?
B=1 + 0.1*cos(f*dT + s1); %voltage of bus
If you take T, it is a vector, and where veer B is used, it will be a vector. Check the formula for B.

Sign in to comment.

Answers (1)

James Tursa
James Tursa on 27 Jun 2022
Edited: James Tursa on 27 Jun 2022
Why are you replacing the sin(y-c) term with approximations? Why haven't you just programmed it as is? It appears you have done this:
sin(y-c) = sin(y)*cos(c) - cos(y)*sin(c)
and then for some reason replaced the sin(y) and cos(y) with truncated Taylor Series approximations, e.g.
sin(y) ~ y - y^3/6
cos(y) ~ 1 - y^2/2
Well, if you know that y is going to remain very small then maybe you can get away with this, but I frankly don't see the point of doing it in the first place unless this eventually will be hosted on a system with no trig functions available.
And for B, since it depends on time, you need to calculate it for the specific time involved in a calculation and make it a part of the inputs to the function handles that depend on it inside the loop. E.g.,
G=@(X1,X2,T,B) etc. % <-- Added B to the input argument list
Although your equations as written don't appear to explicitly depend on T, you are not adjusting T in your RK4 scheme and B depends on this.
Finally, you have two different G's in your code. One is a double variable and the other is a function handle. I suppose this will work since the G in the function handle body gets used prior to the function handle being assigned the same name. But it is confusing code at best. I would advise changing the name of one of those G's. They are both annotated with the same comment "voltage of bus". I suspect that "1" in the B formula is the G=1 from above? If so, use the (changed) variable name rather than a "magic" number. Similar comment for the 0.1 in the B formula.
An outline of what needs to happen with your code is:
G = @(X1,X2,T,B) -(a.*D)./(2.*H1).*X2 + (a.*P)./(2.*H1) -(G.*B.*a)./(2.*H1.*z).*((X1 - (((X1).^3)./6)).*cos(c) -(1-(((X1).^2)./2).*sin(c))); % added B input
B = @(T) 1 + 0.1*cos(f*T + s1); % voltage of bus a function of time
:
for i=1:N-1
K1=F(X2(i));
L1=G(X1(i),X2(i),T(i),B(T(i))); % altered for time
K2=F(X2(i)+((H/2)*L1));
L2=G(X1(i)+((H/2)*K1),X2(i)+((H/2)*L1),T(i)+H/2,B(T(i)+H/2)); % altered for time
K3=F(X2(i)+((H/2)*L2));
L3=G(X1(i)+((H/2)*K2),X2(i)+((H/2)*L2),T(i)+H/2,B(T(i)+H/2)); % altered for time
K4=F(X2(i)+((H)*L3));
L4=G(X1(i)+((H)*K3),X2(i)+((H)*L3),T(i)+H,B(T(i)+H)); % altered for time
X1(i+1)=X1(i)+(H/6)*(K1+2*K2+2*K3+K4);
X2(i+1)=X2(i)+(H/6)*(L1+2*L2+2*L3+L4);
end
You could cut down on some of this code by using a state vector instead of individual variables.

Categories

Find more on Measurements and Testbenches 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!