Error using sym/subs>normalize (line 226) Inconsistency between sizes of second and third arguments.

49 vues (au cours des 30 derniers jours)
Hi, I have an error when calculating equations for Newton method (middle end of code). That Inconsistency between sizes of second and third arguments is refered to subs?
How to deal with that? xout and lambda theta have the same size
Error using sym/subs>normalize (line 226)
Inconsistency between sizes of second and third arguments.
Error in sym/subs>mupadsubs (line 157)
[X2,Y2,symX,symY] = normalize(X,Y); %#ok
Error in sym/subs (line 145)
G = mupadsubs(F,X,Y);
Error in LLC_PGP>Jacobian2D (line 180)
g = subs(f, {lambda theta}, xout); %g(x)
Error in LLC_PGP (line 163)
[n, xout, tol] = Jacobian2D(z,J,x0,tol);
My code below:
% LLC DESIGNER BASED ON OPTIMAL PEAK GAIN PLACEMENT METHOD
%SET PARAMETERS
Vimin=280; %minimum input voltage, where peak gain is placed
Vout=12; %output voltage
RL = 0.24; %max load
Np=16; %primary transformer turns number
Ns=1; %secondary transformer turns number
N=Np./Ns;
Cr=12e-9; %will be calculated later
fsmin=100e3; % Minimum switching frequency from selected transformer Bsat
wsmin=2*pi*fsmin;
%wr=1/(sqrt(Lr*Cr)); %resonant frequency 1 in P mode
%wp=1/(sqrt((Lr+Lp)*Cr)); %resonant frequency 2
t=0;
VCrrating = 2000; %voltag rating of resonant capacitor
%IV. DESIGN ALGORITHM DERIVATION (PN MODE)
%VCr=fVCr(t);
%VCr = k1*cos(wr*t) + k2*sin(wr*t) + Vimin - N*Vo; %eq 10
%iLr = Cr*wr*(-k1*sin(wr*t) + k2*cos(wr*t)); %eq 11
%VCr = k5*cos(wr*t) + k6*sin(wr*t) + Vimin + N*Vo; %eq 12
%iLr = Cr*wr*(-k5*sin(wr*t) + k6*cos(wr*t)); %eq 13
%function VCr = fVCr(t)
% VCr = k1*cos(wr*t) + k2*sin(wr*t) + Vimin - N*Vo;
%end
% 2. Initial conditions: At tf , the resonant current iLr is zero,
% and the Cr voltage is expressed by (6). To solve for the initial condition, t is substituted by 0 in (10) and (11). Then, the
% coefficients k1 and k2 can be solved as
%k1 = (2*N*Cr*RL*Vout*Vimin*fsmin - Cr*RL*(Vimin^2)*fsmin - Vout^2) / (2*RL*fsmin*Cr*Vimin); %eq 14
%k2 = 0; %eq 15
% At tF , the resonant current is also zero, and the Cr voltage
% is expressed by combining (2) and (6). To solve for the initial condition, t is substituted by 0 in (12) and (13). Then, the
% coefficients k5 and k6 can be solved as
%k5 = (2*N*Cr*RL*Vout*Vimin*fsmin + Cr*RL*(Vimin^2)*fsmin - Vout^2) / (-2*RL*fsmin*Cr*Vimin); %eq 16
%k6 = 0; %eq 17
%3) θ and λ Values: According to the definitions in Table II, θ
% and λ have relations with time length
%theta = wr * tfx; %eq 18
%lambda = wr * txf; %eq 19
% Then, the circuit state at tX can be derived from the initial
% condition at tf by substituting (14), (15), and ωr t = q into (10) and (11)
Crmin = (Vout^2 ) / (RL*fsmin*(2*VCrrating-Vimin)*Vimin); %eq 39
theta = acos( (Vimin*(4*(N^2)*Cr*RL*Vout*fsmin - 2*N*Cr*RL*Vimin*fsmin + Vout)) / (2*N*(2*N*Cr*RL*Vout*Vimin*fsmin - Cr*RL*(Vimin^2)*fsmin -(Vout^2))) ); % eq 26
lambda = asin(((2*N*Cr*RL*Vout*Vimin*fsmin - Cr*RL*(Vimin^2)*fsmin -(Vout^2))*sin(theta)) / (2*N*Cr*RL*Vout*Vimin*fsmin + Cr*RL*(Vimin^2)*fsmin -(Vout^2))); % eq 27
% Cr is treated as a known value in (26) and (27). The initial Cr value is
% selected. Based on component stress and then increases with an arbitrary
% interval after each iteration. Each Cr value results in a set of Lr
% and Lp values which provide critically designed peak gain.
% The iteration continues until the boundary between PN and PON modes is reached.
%4) K Value: Lp is clamped by the output capacitor in both
%P and N intervals; therefore, its current slope is linear. The Lp
%current at tf and tF are symmetrical. Also, the Lp current and
%Lr current are equal at tX . Above relations can be described as
% follows:
K = (-N*Cr*RL*Vout*Vimin*fsmin*(lambda + theta)) / ((2*N*Cr*RL*Vout*Vimin*fsmin - Cr*RL*(Vimin^2)*fsmin -(Vout^2))*sin(theta)); % eq 31
% 5) Inductor Values:
wr = 2*lambda*fsmin + 2*theta*fsmin; % eq 33
Lr = 1/(Cr*wr^2); % eq 34
Lp = K*Lr; % eq 35
%6) Checking vCr,tX Margin: When the margin "VCrtXmargin" is less than 0, the design result is invalid, and the design search should continue in PON mode.
VCrtXmargin = (((2*N*Cr*RL*Vout*Vimin*fsmin - Cr*RL*(Vimin^2)*fsmin -(Vout^2))*cos(theta)) / (2*RL*fsmin*Cr*Vimin)) - N*Vout - ((N*Vout*(K+1))/K); %eq 38
% V. DESIGN ALGORITHM DERIVATION (PON MODE)
% 1) Use the last Cr value that has been failed in the PN mode
% as the starting point of the PON mode design search. The
% θ and λ values of the last successful PN mode solution
% are used as the first set of initial values to solve the PON
% mode equations.
%2-d Newtonian's Method for theta and lambda calculation
lambda_init = lambda;
theta_init = theta;
syms lambda theta;
z1 = -2*N*Cr*RL*Vout*Vimin*fsmin*cos(lambda)*lambda ...
- 2*N*Cr*RL*Vout*Vimin*fsmin*cos(lambda)*theta ...
+ 2*N*sin(theta)*Cr*RL*Vout*Vimin*fsmin ...
+ 2*N*sin(lambda)*Cr*RL*Vout*Vimin*fsmin ...
- Cr*RL*(Vimin^2)*fsmin*cos(lambda)*lambda ...
- Cr*RL*(Vimin^2)*fsmin*cos(lambda)*theta ...
- sin(lambda)*(Vout^2) - sin(theta)*Cr*RL*(Vimin^2)*fsmin ...
+ sin(lambda)*Cr*RL*(Vimin^2)*fsmin ...
+ (Vout^2)*cos(lambda)*lambda + (Vout^2)*cos(lambda)*theta - sin(theta)*(Vout^2);
%eq 56
z2 = 4*(2*N*sin(theta)* lambda *Cr*RL* Vout*Vimin*fsmin ...
+ 2*N*sin(theta) *theta*Cr*RL* Vout*Vimin*fsmin ...
- 2*N*sin(lambda) *lambda*Cr*RL* Vout*Vimin*fsmin ...
- 2*N*sin(lambda) *theta*Cr*RL* Vout*Vimin*fsmin ...
+ 4*N*cos(theta) *Cr*RL* Vout*Vimin*fsmin ...
- 4*N*Cr*RL* Vout*Vimin*fsmin*cos(lambda) ...
- sin(theta)*lambda*Cr*RL*(Vimin^2)*fsmin ...
- sin(theta)*theta*Cr*RL*(Vimin^2)*fsmin ...
- sin(lambda)*lambda*Cr*RL*(Vimin^2)*fsmin ...
- sin(lambda)*theta*Cr*RL*(Vimin^2)*fsmin ...
- 2*cos(theta)*Cr*RL*(Vimin^2)*fsmin ...
- 2*Cr*RL*(Vimin^2)*fsmin*cos(lambda) ...
+ 4*Cr*RL*(Vimin^2)*fsmin - sin(theta)*lambda*(Vout^2) ...
- sin(theta)*theta*(Vout^2) + sin(lambda)*lambda*(Vout^2) ...
+ sin(lambda)*theta*(Vout^2) - 2*cos(theta)*(Vout^2) ...
+ 2*(Vout^2)*cos(lambda))*N*Vout*RL*fsmin*(lambda+theta) ...
- 4*(2*N*Cr*RL* Vout*Vimin*fsmin*cos(lambda)* lambda ...
+ 2*N*Cr*RL* Vout*Vimin*fsmin*cos(lambda)* theta ...
- 2*N*sin(theta)*Cr*RL*Vout*Vimin*fsmin ...
- 2*N*sin(lambda)*Cr*RL*Vout*Vimin*fsmin ...
+ Cr*RL*(Vimin^2)*fsmin*cos(lambda)*lambda ...
+ Cr*RL*(Vimin^2)*fsmin*cos(lambda)*theta ...
+ sin(theta)*Cr*RL*(Vimin^2)*fsmin ...
- sin(lambda)*Cr*RL*(Vimin^2)*fsmin ...
- (Vout^2)*cos(lambda)*lambda - (Vout^2)*cos(lambda)*theta + sin(theta)*(Vout^2) ...
+ sin(lambda)*(Vout^2) + 2*lambda*(Vout^2) + 2*theta*(Vout^2))*RL*fsmin*Vimin;
%eq 60
z = [z1; z2];
% Jacobian matrix
J = [diff(z1, 'lambda') diff(z1, 'theta'); diff(z2, 'lambda') diff(z2, 'theta')];
%initial vector
x0 = [lambda_init; theta_init];
% Tolerance:
tol = 1e-10;
%Calculate
[n, xout, tol] = Jacobian2D(z,J,x0,tol);
disp('n = '); disp(n); disp('lambda,theta = '); disp(xout); disp('tol = '); disp(tol);
function [n, xout, tol] = Jacobian2D(f, J, x0, tol)
syms lambda theta;
ep = 1;
n = 0;
xout = x0';
%Computation loop
while ep > tol && n <100
g = subs(f, {lambda theta}, xout); %g(x)
ep = abs(g(1)) + abs(g(2)); %error
Jg = subs(J, {lambda theta}, xout); %Jacobian
yout = xout - Jg\g; %iterate
xout = yout; %update x
n = n + 1; %counter+1
end
xout = double (xout');

Réponse acceptée

Walter Roberson
Walter Roberson le 25 Mai 2022
Jg = subs(J, {lambda theta}, xout); %Jacobian
When you use {} around the list of values to replace, then the third parameter must be a cell array the same size.
You should have coded
Jg = subs(J, [lambda theta], xout); %Jacobian
  3 commentaires
Walter Roberson
Walter Roberson le 26 Mai 2022
You need more iterations to get the desired error level.
% LLC DESIGNER BASED ON OPTIMAL PEAK GAIN PLACEMENT METHOD
%SET PARAMETERS
Vimin=280; %minimum input voltage, where peak gain is placed
Vout=12; %output voltage
RL = 0.24; %max load
Np=16; %primary transformer turns number
Ns=1; %secondary transformer turns number
N=Np./Ns;
Cr=12e-9; %will be calculated later
fsmin=100e3; % Minimum switching frequency from selected transformer Bsat
wsmin=2*pi*fsmin;
%wr=1/(sqrt(Lr*Cr)); %resonant frequency 1 in P mode
%wp=1/(sqrt((Lr+Lp)*Cr)); %resonant frequency 2
t=0;
VCrrating = 2000; %voltag rating of resonant capacitor
%IV. DESIGN ALGORITHM DERIVATION (PN MODE)
%VCr=fVCr(t);
%VCr = k1*cos(wr*t) + k2*sin(wr*t) + Vimin - N*Vo; %eq 10
%iLr = Cr*wr*(-k1*sin(wr*t) + k2*cos(wr*t)); %eq 11
%VCr = k5*cos(wr*t) + k6*sin(wr*t) + Vimin + N*Vo; %eq 12
%iLr = Cr*wr*(-k5*sin(wr*t) + k6*cos(wr*t)); %eq 13
%function VCr = fVCr(t)
% VCr = k1*cos(wr*t) + k2*sin(wr*t) + Vimin - N*Vo;
%end
% 2. Initial conditions: At tf , the resonant current iLr is zero,
% and the Cr voltage is expressed by (6). To solve for the initial condition, t is substituted by 0 in (10) and (11). Then, the
% coefficients k1 and k2 can be solved as
%k1 = (2*N*Cr*RL*Vout*Vimin*fsmin - Cr*RL*(Vimin^2)*fsmin - Vout^2) / (2*RL*fsmin*Cr*Vimin); %eq 14
%k2 = 0; %eq 15
% At tF , the resonant current is also zero, and the Cr voltage
% is expressed by combining (2) and (6). To solve for the initial condition, t is substituted by 0 in (12) and (13). Then, the
% coefficients k5 and k6 can be solved as
%k5 = (2*N*Cr*RL*Vout*Vimin*fsmin + Cr*RL*(Vimin^2)*fsmin - Vout^2) / (-2*RL*fsmin*Cr*Vimin); %eq 16
%k6 = 0; %eq 17
%3) θ and λ Values: According to the definitions in Table II, θ
% and λ have relations with time length
%theta = wr * tfx; %eq 18
%lambda = wr * txf; %eq 19
% Then, the circuit state at tX can be derived from the initial
% condition at tf by substituting (14), (15), and ωr t = q into (10) and (11)
Crmin = (Vout^2 ) / (RL*fsmin*(2*VCrrating-Vimin)*Vimin); %eq 39
theta = acos( (Vimin*(4*(N^2)*Cr*RL*Vout*fsmin - 2*N*Cr*RL*Vimin*fsmin + Vout)) / (2*N*(2*N*Cr*RL*Vout*Vimin*fsmin - Cr*RL*(Vimin^2)*fsmin -(Vout^2))) ); % eq 26
lambda = asin(((2*N*Cr*RL*Vout*Vimin*fsmin - Cr*RL*(Vimin^2)*fsmin -(Vout^2))*sin(theta)) / (2*N*Cr*RL*Vout*Vimin*fsmin + Cr*RL*(Vimin^2)*fsmin -(Vout^2))); % eq 27
% Cr is treated as a known value in (26) and (27). The initial Cr value is
% selected. Based on component stress and then increases with an arbitrary
% interval after each iteration. Each Cr value results in a set of Lr
% and Lp values which provide critically designed peak gain.
% The iteration continues until the boundary between PN and PON modes is reached.
%4) K Value: Lp is clamped by the output capacitor in both
%P and N intervals; therefore, its current slope is linear. The Lp
%current at tf and tF are symmetrical. Also, the Lp current and
%Lr current are equal at tX . Above relations can be described as
% follows:
K = (-N*Cr*RL*Vout*Vimin*fsmin*(lambda + theta)) / ((2*N*Cr*RL*Vout*Vimin*fsmin - Cr*RL*(Vimin^2)*fsmin -(Vout^2))*sin(theta)); % eq 31
% 5) Inductor Values:
wr = 2*lambda*fsmin + 2*theta*fsmin; % eq 33
Lr = 1/(Cr*wr^2); % eq 34
Lp = K*Lr; % eq 35
%6) Checking vCr,tX Margin: When the margin "VCrtXmargin" is less than 0, the design result is invalid, and the design search should continue in PON mode.
VCrtXmargin = (((2*N*Cr*RL*Vout*Vimin*fsmin - Cr*RL*(Vimin^2)*fsmin -(Vout^2))*cos(theta)) / (2*RL*fsmin*Cr*Vimin)) - N*Vout - ((N*Vout*(K+1))/K); %eq 38
% V. DESIGN ALGORITHM DERIVATION (PON MODE)
% 1) Use the last Cr value that has been failed in the PN mode
% as the starting point of the PON mode design search. The
% θ and λ values of the last successful PN mode solution
% are used as the first set of initial values to solve the PON
% mode equations.
%2-d Newtonian's Method for theta and lambda calculation
lambda_init = lambda;
theta_init = theta;
syms lambda theta;
z1 = -2*N*Cr*RL*Vout*Vimin*fsmin*cos(lambda)*lambda ...
- 2*N*Cr*RL*Vout*Vimin*fsmin*cos(lambda)*theta ...
+ 2*N*sin(theta)*Cr*RL*Vout*Vimin*fsmin ...
+ 2*N*sin(lambda)*Cr*RL*Vout*Vimin*fsmin ...
- Cr*RL*(Vimin^2)*fsmin*cos(lambda)*lambda ...
- Cr*RL*(Vimin^2)*fsmin*cos(lambda)*theta ...
- sin(lambda)*(Vout^2) - sin(theta)*Cr*RL*(Vimin^2)*fsmin ...
+ sin(lambda)*Cr*RL*(Vimin^2)*fsmin ...
+ (Vout^2)*cos(lambda)*lambda + (Vout^2)*cos(lambda)*theta - sin(theta)*(Vout^2);
%eq 56
z2 = 4*(2*N*sin(theta)* lambda *Cr*RL* Vout*Vimin*fsmin ...
+ 2*N*sin(theta) *theta*Cr*RL* Vout*Vimin*fsmin ...
- 2*N*sin(lambda) *lambda*Cr*RL* Vout*Vimin*fsmin ...
- 2*N*sin(lambda) *theta*Cr*RL* Vout*Vimin*fsmin ...
+ 4*N*cos(theta) *Cr*RL* Vout*Vimin*fsmin ...
- 4*N*Cr*RL* Vout*Vimin*fsmin*cos(lambda) ...
- sin(theta)*lambda*Cr*RL*(Vimin^2)*fsmin ...
- sin(theta)*theta*Cr*RL*(Vimin^2)*fsmin ...
- sin(lambda)*lambda*Cr*RL*(Vimin^2)*fsmin ...
- sin(lambda)*theta*Cr*RL*(Vimin^2)*fsmin ...
- 2*cos(theta)*Cr*RL*(Vimin^2)*fsmin ...
- 2*Cr*RL*(Vimin^2)*fsmin*cos(lambda) ...
+ 4*Cr*RL*(Vimin^2)*fsmin - sin(theta)*lambda*(Vout^2) ...
- sin(theta)*theta*(Vout^2) + sin(lambda)*lambda*(Vout^2) ...
+ sin(lambda)*theta*(Vout^2) - 2*cos(theta)*(Vout^2) ...
+ 2*(Vout^2)*cos(lambda))*N*Vout*RL*fsmin*(lambda+theta) ...
- 4*(2*N*Cr*RL* Vout*Vimin*fsmin*cos(lambda)* lambda ...
+ 2*N*Cr*RL* Vout*Vimin*fsmin*cos(lambda)* theta ...
- 2*N*sin(theta)*Cr*RL*Vout*Vimin*fsmin ...
- 2*N*sin(lambda)*Cr*RL*Vout*Vimin*fsmin ...
+ Cr*RL*(Vimin^2)*fsmin*cos(lambda)*lambda ...
+ Cr*RL*(Vimin^2)*fsmin*cos(lambda)*theta ...
+ sin(theta)*Cr*RL*(Vimin^2)*fsmin ...
- sin(lambda)*Cr*RL*(Vimin^2)*fsmin ...
- (Vout^2)*cos(lambda)*lambda - (Vout^2)*cos(lambda)*theta + sin(theta)*(Vout^2) ...
+ sin(lambda)*(Vout^2) + 2*lambda*(Vout^2) + 2*theta*(Vout^2))*RL*fsmin*Vimin;
%eq 60
z = [z1; z2];
% Jacobian matrix
J = [diff(z1, 'lambda') diff(z1, 'theta'); diff(z2, 'lambda') diff(z2, 'theta')];
%initial vector
x0 = [lambda_init; theta_init];
% Tolerance:
tol = 1e-10;
%Calculate
[n, xout, tol, ep] = Jacobian2D(z,J,x0,tol);
disp('n = '); disp(n); disp('lambda,theta = '); disp(xout); disp('tol = '); disp(tol);
n = 100 lambda,theta = 1.0818 2.5575 tol = 1.0000e-10
fprintf('Err = %g\n', double(ep));
Err = 1.27516e-06
function [n, xout, tol, ep] = Jacobian2D(f, J, x0, tol)
syms lambda theta;
ep = 1;
n = 0;
xout = x0;
%Computation loop
while ep > tol && n <100
%whos
g = subs(f, [lambda; theta], xout); %g(x)
ep = abs(g(1)) + abs(g(2)); %error
Jg = subs(J, [lambda; theta], xout); %Jacobian
yout = xout - double(Jg)\double(g); %iterate
xout = yout; %update x
n = n + 1; %counter+1
end
xout = double (xout');
end
Marcin Fisior
Marcin Fisior le 26 Mai 2022
Great, thanks. So you also removed the transposition of x0, and added the double conversion of Jg and g.
Now its working.

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Numeric Types dans Help Center et File Exchange

Tags

Produits


Version

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by