Clear Filters
Clear Filters

Solve systems of 6 non-linear equations

14 views (last 30 days)
I am trying to solve a sytem of 6 non-linear equations. I used vpasolve. One solution it gave me is I1=I2, V1=V2, and hence, my deltaT2 is roughly 0. So, I set the starting values of all to 1e-10 so that my deltaT2 will be non-zero. But after around one minute, I got an empty solution. I believe there should be a solution set of values for this system. So I wonder what could go wrong with my code/solution.
To make those equations I am solving more sense, I attached the screenshot I got from this paper below. Equation 2-3 are derived by using two points at t0 and t1 lying on one circle. Equation 4-5 are derived by using two points at t1 and t2 lying on an eclipse.
The only solution I got from my code so far is that V1=V2 and I1=I2. But as you can see, I want V1 and V2 to have different values. I1 and I2 should have differenet values as well.
clear;clc;close all;
%% Parameters
Vin = 400; Lr = 55e-6; Cr = 24e-9; Lm = 280e-6;
Zo = sqrt(Lr/Cr); Vo = 12; n = Vin/2/Vo;
wo = 1/sqrt(Lr*Cr); fo = wo/(2*pi);
Iheavy = 15; Rheavy = Vo/Iheavy;
Ilight = 5; Rlight = Vo/Ilight;
fs = 120e3; Ts = 1/fs;
Vin = 350;
Zr1 = sqrt(Lr/Cr); wr1 = sqrt(1/(Lr*Cr));
Zr2 = sqrt((Lr+Lm)/Cr); wr2 = sqrt(1/((Lr+Lm)*Cr));
R = Rheavy;
syms I1 V1 V2 I2 deltaT1 deltaT2 real
eq1 = I1 + I2 - (n*Vo/Lm)*deltaT1;
eq2 = (V1 - Vin + n*Vo)^2 + (I1*Zr1)^2 - (V2 - Vin + n*Vo)^2 - (I2*Zr1)^2;
eq3 = deltaT1 - (pi-atan(I2*Zr1/(V2-Vin+n*Vo)) + atan(I1*Zr1/(V1+Vin-n*Vo)))/wr1;
eq4 = (V1 - Vin)^2 + (I1*Zr2)^2 - (V2 - Vin)^2 - (I2*Zr2)^2;
eq5 = deltaT2 - (atan(I1*Zr2/(Vin-V1)) - atan(I2*Zr2/(Vin-V2)))/wr2;
eq6 = Vo*(deltaT1 + deltaT2)/(n*R) - (V1+V2)*Cr - ((I1-I2)/2)*deltaT1;
eqs = [eq1,eq2,eq3,eq4,eq5,eq6];
starting_value = [1e-10 Inf; 1e-10 Inf; 1e-10 Inf; ...
1e-10 Inf; 1e-10 Inf; 1e-10 Inf];
[I1,V1,V2,I2,deltaT1,deltaT2] = vpasolve(eqs,[I1,V1,V2,I2,deltaT1,deltaT2],...
starting_value);
result = double([I1,V1,V2,I2,deltaT1,deltaT2])

Accepted Answer

Tra Yen Nhu Phan
Tra Yen Nhu Phan on 25 Oct 2023
Edited: Tra Yen Nhu Phan on 25 Oct 2023
I figured the problem out. The paper has a small syntax in their system of equations. Therefore, MATLAB cannot give the result that I want. The mistake was at equation 2 that Vin should be -Vin. Also, instead of using my own parameters, I used the paper's parameters, and everything just works perfectly. Thanks all for helping me this out.
Below is the code that works for me. The modification includes:
1) Correct the error in the equation system
2) Use parameters in the paper instead
3) Add the range for all parameters as positive numbers
clear;clc; close all;
Vo = 48; Po = 7200; figure;
N = 18; Lm = 111.4e-6; Lr = 19.18e-6;
Cr = 20.25e-9; Co = 1e-3;
Zr1 = sqrt(Lr/Cr); wr1 = 1/sqrt(Lr*Cr);
fr1 = wr1/(2*pi); Tr1 = 1/fr1;
Zr2 = sqrt((Lr + Lm)/Cr); wr2 = sqrt(1/((Lr + Lm)*Cr));
fr2 = wr2/(2*pi); Tr2 = 1/fr2;
select = 1;
Q1 = 0.1831; Q2 = 0.3662;
if (select == 1)
Q = Q1;
else
Q = Q2;
end
Rac = Zr1/Q; Ro = Rac*(pi^2)/8/(N^2);
%% Below Resonance
syms I1 V1 I2 V2 deltaT1 deltaT2 fs
Vin = 800;
eq1 = (I1 + I2) -(N*Vo/Lm)*deltaT1;
eq2 = (-V1 - Vin + N*Vo)^2 +(I1*Zr1)^2 - (V2 - Vin + N*Vo)^2 - (I2*Zr1)^2 ;
eq3 = deltaT1 - ((pi - atan(I2*Zr1/(V2-Vin+N*Vo)) + atan(I1*Zr1/(V1+Vin-N*Vo)))/wr1);
eq4 = (V1 - Vin)^2 +(I1*Zr2)^2 - (V2 - Vin)^2 - (I2*Zr2)^2 ;
eq5 = deltaT2 - ((atan(I1*Zr2/(Vin - V1)) - atan(I2*Zr2/(Vin - V2)))/wr2);
eq6 = (Vo*(deltaT1 + deltaT2)/(N*Ro)) - (V1+V2)*Cr - (deltaT1*(I1-I2))/2;
eqs = [eq1,eq2,eq3,eq4,eq5,eq6];
range = [0, Inf; % I1
0, Inf; % V1
0, Inf; % I2
0, Inf; % V2
0, Inf; % deltaT1
0, Inf]; % deltaT2
[I1,V1,V2,I2,deltaT1,deltaT2] = vpasolve(eqs,[I1,V1,V2,I2,deltaT1,deltaT2],range)
I1 = 
8.483487211088991203528612806318
V1 = 
253.64132442119089964849356851349
V2 = 
163.29575706350444322315252328323
I2 = 
7.4427274570926256928348488481962
deltaT1 = 
0.0000020534494375410093334689879526685
deltaT2 = 
0.00000022936567635361382499350827474835
fs = 1/(deltaT1*2 + deltaT2*2)
fs = 
219027.81217659331611383356783722
  5 Comments
Tra Yen Nhu Phan
Tra Yen Nhu Phan on 25 Oct 2023
I absolutely can do that. Thanks for letting me know what I should do to express my appreciation for all's help.

Sign in to comment.

More Answers (0)

Categories

Find more on Interactive Control and Callbacks in Help Center and File Exchange

Products


Release

R2019a

Community Treasure Hunt

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

Start Hunting!