Nonlinear Inductor Matlab Code

In a CP-line model, consider the case in which the transmission line is feeding a nonlinear inductor:

• Assume a simple two slope model for the magnetization characteristic of the inductor.

• Use only the trapezoidal rule of integration (you may want to consider CDA if you encounter numerical oscillations).

• Assume zero initial conditions everywhere.

Use Matlab code for this question, and no simulink just code.

Also below is matlab code for the cp-line model without the nonlinear inductor

DT = 0.000037;
Tmax=0.02;
Topen=0.008;

Zc=288.675;
R1=3;
L1=350/1e3;
C1=10/1e9;
C=12/1e9;
R=0.5;
L=1/1e3;
A=230*1e3;
w=2*pi*60;
tau=0.000001732;

Kmax=round(Tmax/DT);
Kopen=ceil(Topen/DT)-1; % begin checking current direction at Kopen

v1(1)=0; % initial condition
v2(1)=0;
v3(1)=0;
v4(1)=0;
v5(1)=0;
v6(1)=0;
v7(1)=0;
i12(1)=0;
i70(1)=0;
vsw(1)=0;
isw(1)=0;
h12(1)=0;
h20(1)=0;
h40(1)=0;
h50(1)=0;
h60(1)=0;
h70(1)=0;

gs12=1/(R1+2*L1/DT);
g20=2*C1/DT;
G1=[gs12+g20+4/R   -4/R       0         0         0;
    -4/R           4/R+1/Zc   0         0         0;
    0              0          2/R+1/Zc  -2/R      0;
    0              0          -2/R      2/R+1/Zc  0;
    0              0          0         0         4/R+1/Zc];
G2=[gs12+g20   0      0         0          0          0;
    0          4/R    -4/R      0          0          0;
    0          -4/R   4/R+1/Zc  0          0          0;
    0          0      0         2/R+1/Zc   -2/R       0;
    0          0      0         -2/R       2/R+1/Zc   0;
    0          0      0         0          0          4/R+1/Zc];

i=2;
flag = 1; % flag for closed switch. if 0, then switch is opened
%switch closed
while (i<=Kmax+1 & flag) % history terms
    h12(i)=(-1*(v1(i-1)-v2(i-1))-i12(i-1)*(2*L1/DT-R1))*gs12;
    h20(i)=v2(i-1)*2*g20-h20(i-1);
    if (i-1)*DT <= tau
        h40(i)=0;
        h50(i)=0;
        h60(i)=0;
        h70(i)=0;
    else
        upper=ceil(i-tau/DT);
        lower=floor(i-tau/DT);
        h40(i)=((i-tau/DT-lower)*(v5(upper)-v5(lower))+v5(lower))*2/Zc-((i-tau/DT-lower)*(h50(upper)-h50(lower))+h50(lower));
        h50(i)=((i-tau/DT-lower)*(v4(upper)-v4(lower))+v4(lower))*2/Zc-((i-tau/DT-lower)*(h40(upper)-h40(lower))+h40(lower));
        h60(i)=((i-tau/DT-lower)*(v7(upper)-v7(lower))+v7(lower))*2/Zc-((i-tau/DT-lower)*(h70(upper)-h70(lower))+h70(lower));
        h70(i)=((i-tau/DT-lower)*(v6(upper)-v6(lower))+v6(lower))*2/Zc-((i-tau/DT-lower)*(h60(upper)-h60(lower))+h60(lower));
    end
    v1(i)=A*cos(w*(i-1)*DT);
    Hist1=[-1*h12(i)+h20(i); h40(i); h50(i); h60(i); h70(i)];
    Hist2=[-1*gs12; 0; 0; 0; 0];
    X=G1\(Hist1-Hist2*v1(i));
    v2(i)=X(1);
    v4(i)=X(2);
    v5(i)=X(3);
    v6(i)=X(4);
    v7(i)=X(5);
    v3(i)=v2(i);
    i12(i)=(v1(i)-v2(i))*gs12-h12(i);
    i70(i)=v7(i)*4/R;
    vsw(i)=v2(i)-v3(i);
    isw(i)=h20(i)-v2(i)*g20+(-h12(i)-(v2(i)-v1(i))*gs12);
    if i>Kopen+1
        if sign(isw(i))~=sign(isw(i-1)) %if the current direction changes, open the switch
            flag=0;
        end
    end
    i=i+1;
end
K=i;

%switch open
while (i<=Kmax+1 & ~flag) % history terms
    h12(i) = (-1*(v1(i-1)-v2(i-1))-i12(i-1)*(2*L1/DT-R1))*gs12;
    h20(i) = v2(i-1)*2*g20-h20(i-1);
    if(i-1)*DT < tau
        h40(i)=0;
        h50(i)=0;
        h60(i)=0;
        h70(i)=0;
    else
        upper = ceil(i-tau/DT);
        lower = floor(i-tau/DT);
        h40(i)=((i-tau/DT-lower)*(v5(upper)-v5(lower))+v5(lower))*2/Zc-((i-tau/DT-lower)*(h50(upper)-h50(lower))+h50(lower));
        h50(i)=((i-tau/DT-lower)*(v4(upper)-v4(lower))+v4(lower))*2/Zc-((i-tau/DT-lower)*(h40(upper)-h40(lower))+h40(lower));
        h60(i)=((i-tau/DT-lower)*(v7(upper)-v7(lower))+v7(lower))*2/Zc-((i-tau/DT-lower)*(h70(upper)-h70(lower))+h70(lower));
        h70(i)=((i-tau/DT-lower)*(v6(upper)-v6(lower))+v6(lower))*2/Zc-((i-tau/DT-lower)*(h60(upper)-h60(lower))+h60(lower));
    end
    v1(i)=A*cos(w*(i-1)*DT);
    Hist1=[-1*h12(i)+h20(i); 0; h40(i); h50(i); h60(i); h70(i)];
    Hist2=[-1*gs12; 0; 0; 0; 0; 0];
    X=G2\(Hist1-Hist2*v1(i));
    v2(i) = X(1);
    v3(i) = X(2);
    v4(i) = X(3);
    v5(i) = X(4);
    v6(i) = X(5);
    v7(i) = X(6);
    
    i12(i)=(v1(i)-v2(i))*gs12-h12(i);
    i70(i)=v7(i)*4/R;
    vsw(i)=v2(i)-v3(i);
    isw(i)=0;
    i=i+1;
end


x1 = linspace(0,.02,length(v1));
plot(x1,v1)
xlabel('Time(s)');
ylabel('V1 (volts)');
title('V1');

figure;

x2 = linspace(0,.02,length(v2));
plot(x2,v2)
xlabel('Time(s)');
ylabel('V2 (volts)');
title('V2');

figure;

x3 = linspace(0,.02,length(vsw));
plot(x3,vsw)
xlabel('Time(s)');
ylabel('VSW (volts)');
title('Voltage across the switch');

figure;

x4 = linspace(0,.02,length(i70));
plot(x4,i70)
xlabel('Time(s)');
ylabel('I70 (amps)');
title('Short Circuit Current');