¿Cómo resolver un circuito boost (o cualquier circuito RLC conmutado) en el dominio del tiempo?

Question

Necesito modelar un circuito boost con la siguiente topología:

simular este circuito - Esquema creado con CircuitLab

Me interesa la tensión de la resistencia, así que encontré la ecuación diferencial de cada estado y las resolví por el método implícito de Euler, alternando la ecuación a resolver según el cambio de estado.

He comparado los resultados con un simulador de circuitos. Cada ecuación diferencial es correcta por sí sola, pero cuando las pongo juntas, una tras otra, formando la dinámica conmutada de boost, el resultado es erróneo.

Estos son mis pasos:

Ecuación diferencial del "estado ON" (sólo la descarga del condensador)

$$\frac{dV_{out}}{dt} = -\frac{V_{out}}{RC}$$

Ecuación diferencial "estado OFF

$$\frac{d^2V_{out}}{dt^2} + \frac{1}{RC}\frac{dV_{out}}{dt} + \frac{1}{LC}V_{out} = \frac{V_{in}}{LC}$$

Este es mi código (MATLAB):

R = 10; 
L = 1e-5; 
C = 1e-4;
Vin = 10; 

h = 1e-6;   % time step
sizeVector = 10000; % number of steps
t = 0:h:sizeVector*h-h;  % time vector

Vout = zeros(1,sizeVector); % output voltage vector
dVout = zeros(1,sizeVector); % first derivative vector

freqSw = 5e3; % switching frequency
D = zeros(1,sizeVector); % State of switch for each simulation step (duty cycle = 0.5)
D = sign(cos(2*pi*freqSw.*t + pi/2 + 0.001));
D(D<0) = 0;

for k=2:sizeVector
    if D(k) == 0
        % Off state
        Vout(k) = (Vout(k-1) + h*((dVout(k-1) + (h/(L*C))*Vin)/(1 + h/(R*C))))/(1 + (h^2/(L*C))/(1 + h/(R*C)));
        dVout(k) = (dVout(k-1) + h*(Vin/(L*C) - Vout(k)/(L*C)))/(1 + h/(R*C));
    end
    if D(k) == 1
        % On state
        Vout(k) = Vout(k-1)/(1+h/(R*C));
        dVout(k) = (Vout(k)-Vout(k-1))/h; % just for convenience
    end

end

figure(1)
plot(t,Vout)
xlabel('Time (s)')
ylabel('Voltage (V)')
hold on

Estas son las comparaciones entre los resultados de PSIM (rojo) y mi código en MATLAB (azul). (Para la descarga del condensador, he considerado en ambos casos una tensión inicial de 10V)

Estado abierto:

Estado cerrado:

Convertidor Boost como en el código anterior:

Nunca puedo conseguir una tensión de estado estable superior a mi tensión de entrada (10V, en este caso), pero la forma de onda parece coherente con cada uno de los estados de forma de onda por separado..

¿Debe haber un error en mi código o me estoy perdiendo algún punto importante? ¿Es posible este tipo de solución?

Editado después de las respuestas:

Siguiendo las sugerencias, cambié la ecuación diferencial de segundo orden por un sistema de dos ecuaciones de primer orden, para poder calcular la corriente del inductor, lo cual es necesario porque cuando comienza el "estado de apagado", la corriente del inductor debe ser una condición inicial. En el modelo anterior, no estaba considerando el cambio de corriente en el 'estado de encendido'.

El sistema "off state":

$$\begin{equation} \begin{bmatrix} \frac{dV_{out}}{dt}\\ \frac{dI_{ind}}{dt} \end{bmatrix} = \begin{bmatrix} -\frac{1}{RC} & \frac{1}{C}\\ -\frac{1}{L} & 0 \end{bmatrix} | tiempos \begin{bmatrix} V_{out}\\ I_{ind} \end{bmatrix} + \begin{bmatrix} 0\\ \frac{Vin}{L} \end{bmatrix} \end{equation}$$

El sistema 'on state': (en realidad no es un sistema, pero lo mantengo matricial para simplificar)

$$\begin{equation} \begin{bmatrix} \frac{dV_{out}}{dt}\\ \frac{dI_{ind}}{dt} \end{bmatrix} = \begin{bmatrix} -\frac{1}{RC} & 0\\ 0 & 0 \end{bmatrix} | tiempos \begin{bmatrix} V_{out}\\ I_{ind} \end{bmatrix} + \begin{bmatrix} 0\\ \frac{Vin}{L} \end{bmatrix} \end{equation}$$

Todavía hay un problema. De esta manera, no puedo evitar que la corriente del inductor sea negativa en el "estado apagado", por lo que inserté una sentencia if en cada iteración del "estado apagado". Si $I_{ind} < 0$ , resuelve el sistema de nuevo imponiendo $I_{ind} = 0$ . Esto se puede hacer simplemente poniendo a cero los elementos relacionados con $I_{ind}$ en la matriz.

A continuación se muestra la comparación de resultados (PSIM, primero, MATLAB, después) y el código.

R = 10; 
L = 1e-3; 
C = 1e-4;
Vin = 10; % DC input voltage

h = 1e-7;                % step
sizeVector = 120000;     % number of steps
t = 0:h:sizeVector*h-h;  % time vector

X = zeros(2,sizeVector); % state vector 

A = [-1/(R*C) 1/C ; -1/L 0]; % off state matrix
b = [0; Vin/L];              % off state vector

Aaux = [-1/(R*C) 0 ; 0 0];   % off state auxiliar matrix
baux = [0 ; 0];              % off state auxiliar vector

A2 = [-1/(R*C) 0; 0 0];      % on state matrix
b2 = [0; Vin/L];             % on state vector

freqSw = 5e3;               % switching frequency
D = zeros(1,sizeVector); % State of switch for each simulation step (duty cycle = 0.5)
D = sign(cos(2*pi*freqSw.*t + pi/2 + 0.001));
D(D<0) = 0;

for k=2:sizeVector

    if D(k) == 0 % off state 
        X(:,k) = (inv(eye(length(A)) - h*A))*(X(:,k-1) + h*b);
        if X(2,k) < 0 % avoid Iind < 0
            X(:,k) = (inv(eye(length(Aaux)) - h*Aaux))*(X(:,k-1) + h*baux);
        end

    end
    if D(k) == 1 % on state
        X(:,k) = (inv(eye(length(A2)) - h*A2))*(X(:,k-1) + h*b2);
    end
end

figure(1)
plot(t,X(1,:))
hold on
plot(t,X(2,:))
xlabel('Time (s)')
ylabel('Voltage (V)\Current(A)')
legend('V_{out}','I_{in}')
ylim([-6 36])
grid
hold on

Answer 1

Veo al menos una cosa mal en tus ecuaciones (no he revisado tu código):

• Debe haber una ecuación que describa la evolución de la corriente del inductor durante el "estado de encendido".

Probablemente me quedaría con un sistema de dos ecuaciones de primer orden en lugar de intentar escribir una ecuación de segundo orden. No va a tomar Matlab más tiempo para resolverlo de esa manera.

Answer 2

Una forma de analizar la conversión boost es calcular la energía almacenada en el inductor basándose en el tiempo de conmutación, Ton, la corriente, la inductancia I, L cuando el interruptor está cerrado y luego equiparar eso a la energía disipada en la carga en un bucle regulado.

P_in = 0,5 * L I²/T ( Energía/tiempo cuando el interruptor está en ON)
Pout = 0,5 * I * Vin * Toff/T (Power Out)
Suponiendo una conmutación sin pérdidas
tensión de salida $Vo=Vin\sqrt{\dfrac{kR_LT_{on}}{2L}}$ para 0< k< 1 con un bucle regulado