Расчёт переходных процессов методом синтетических схем, страница 4

for k1 = 1:D1               %Обращение всех 0 из матрицы R в ноль inf.

    if R(k1,1)~=0           %Необходимо для поиска минимального

        R1(k1,1)=R(k1,1);   %значения сопротивления.

    else

        R1(k1,1)=inf ;

    end

end

Rmin=min(R1);               % Минимальное значение сопротивления.

L1=zeros(D1,1);

for k2 = 1:D1               %Обращение всех 0 из матрицы L в ноль inf.

    if L(k2,1)~=0           %Необходимо для поиска минимального

        L1(k2,1)=L(k2,1);   %значения индуктивности.

    else

        L1(k2,1)=inf ;

    end

end

Lmin=min(L1);                % Минимальное значение индуктивности.

if (Lmin/Rmax)<(C1max*Rmin)  % Оценка шага интегрирования.

    h=0.1*(Lmin/Rmax);

else

    h=0.1*(C1max*Rmin);

end

if (Lmax/Rmin)>(Cmin*Rmax)   % Нахождение максимального времени расчёта.

    T=5*(Lmax/Rmin);         

else

    T=5*(Cmin*Rmax);

end

V1=NET(:,2);

V2=NET(:,3);

A=zeros(q,p);                % Построение матрицы узловых соединений.

for k3=1:p

    A(V1(k3,1),k3)=1;

    A(V2(k3,1),k3)=-1;

end

A(q,:)=[];

for k4=1:p

    if C(k4,1)~=inf                       % Макромодель RC-цепи.

        g(k4,1)=1/(R(k4,1)+h/(C(k4,1)));  % Проводимость RC-цепи.

        J1(k4,1)=-Nach(k4,1)*g(k4,1);     % Источник тока RC-цепи.

        Ukon(k4,1)=Nach(k4,1);  % Значение напряжений из начальных данных

    elseif L(k4,1)~=0                 % Макромодель RL-цепи.

        g(k4,1)=h/(L(k4,1)+h*R(k4,1));         % Проводимость RL-цепи.

        J2(k4,1)= L(k4,1)/(L(k4,1)+h*R(k4,1)); % Источник тока RL-цепи.

        J1(k4,1)=Nach(k4,1).*J2(k4,1);

        I(k4,1)=Nach(k4,1);     % Значение токов из начальных данных

    elseif R(k4,1)~=0

        g(k4,1)=1/(R(k4,1)); % Проводимость цепи состоящей только из R.

        J1(k4,1)=J(k4,1);

    else

        g(k4,1)=0;          % Проводимость цепи состоящей только из J.

        J1(k4,1)=J(k4,1);

    end

end

H=diag(g);

G=inv(A*H*(A.'));                   %Матрица узловых проводимостей.

Toki=zeros(p,1);

Ugr=zeros(p,1);

M=zeros(1,fix(T/h));

for k5=1:T/h

    B=-A*(J1+H*E);                  %Матрица задающих токов узлов.         

    U0=G*B;                         %Значение потенциалов узлов.       

    Up=(A.')*U0;                    %Внешние напряжения ветвей.    

    U=Up+E;                         %Внутренние напряжения ветвей.    

    I=H*U;                          %Внутренние токи ветвей       

    Ip=I+J1;                        %Внешние токи ветвей       

    Ukon=Up-Ip.*R;                  %Напряжение на конденсаторах и ЕДС.

    for k6=1:p                      % Перерасчёт источников тока на:

        if L(k6,1)~=0               %Цепях с индуктивностью

            J1(k6,1)=Ip(k6,1).*J2(k6,1);

        elseif C(k6,1)~=inf         %Цепях с ёмкостью

            J1(k6,1)=-Ukon(k6,1)*H(k6,k6);

        end

    end

    Toki=[Toki Ip]; %Вывод матрицы токов все ветви и моменты времени.

    Ugr=[Ugr U];    %Вывод матрицы напряжений все ветви и моменты времени.

    M(1,k5)=k5;

end

k=1:1:(T/h+1);    

subplot(1,2,1);

plot(k,Toki(1,:));     %Построение зависимостей токов от времени.

for k10=2:p

    hold on

    plot(k,Toki(k10,:));

end

subplot(1,2,2);       

plot (k,Ugr(1,:));     %Построение зависимостей напряжений от времени.

for k11=2:p

    hold on

    plot (k,Ugr(k11,:));

end

end