% Trabajo de Fin de Grado.
% Grado en Ingeniería Electrónica Industrial.
% Título: Diseño y control de un convertidor Back-to-Back: Análisis teórico
% y validación mediante simulación.
% Autor: José Antonio Ramírez Jurado.
% Tutor: Dr. Manuel Jesús López Sánchez.
% Cádiz, junio 2024.

%--------------------------------------------------------------------------

% Subtítulo: Modelo no lineal. LQR con integral.

% Parámetros constantes del VSC1.
L_1 = 4.1e-3;           % Inductancia (henrios).
R_1 = 284e-3;           % Resistencia (ohmios).
m_t1 = 5;               % Amplitud portadora triangular (voltios).
omega_1 = 377;          % Frecuencia (radianes/segundo).
v_d_1 = 100*sqrt(2);    % Tension pico AC VSC1 (Voltios).
                        % V_pico = V_rms * sqrt(2).

% Parámetros constantes del VSC2.
L_2 = 5.3e-3;           % Inductancia (henrios).
R_2 = 330e-3;           % Resistencia (ohmios).
m_t2 = 5;               % Amplitud portadora triangular (voltios).
omega_2 = 377;          % Frecuencia (radianes/segundo).
v_d_2 = 100*sqrt(2);    % Tension pico AC VSC1 (Voltios).
                        % V_pico = V_rms * sqrt(2).

% Parámetros constantes del sistema.
C_DC = 6.8e-3;          % Capacitancia condensador (Faradios).

% Calculo de la condición de operación del modelo no lineal.
% Condiciones de operación arbitrarias.
V_DC0 = 500; 
i_d0_1=((3000*2)/(3*v_d_1))*1.06522;  
i_q0_1 = 0;
i_d0_2= -i_d0_1/1.06522;
i_q0_2 = 0;

% Variables de control despejadas de las primeras 4 ecuaciones bajo la
% condición de que el sistema está en equilibrio, es decir, las variables
% de estado no varían en el tiempo.
eta_d0_1 = ((i_q0_1*omega_1 + v_d_1/L_1 - (R_1/L_1)*i_d0_1)*(2*m_t1*L_1))/V_DC0
eta_q0_1 = (- i_d0_1*omega_1 - (R_1/L_1)*i_q0_1)*(2*m_t1*L_1)/V_DC0 
eta_d0_2 = ((i_q0_2*omega_2 + v_d_2/L_2 - (R_2/L_2)*i_d0_2)*(2*m_t2*L_2))/V_DC0  
eta_q0_2 = ((- i_d0_2*omega_2 - (R_2/L_2)*i_q0_2)*(2*m_t2*L_2))/V_DC0

% Tras los cálculos anteriores, la variación de V_DC tiene que ser nula o
% prácticamente nula.
dV_DC_d= 3/(4*C_DC)*((i_d0_1*eta_d0_1 + i_q0_1*eta_q0_1)/m_t1) + ...
3/(4*C_DC)*((i_d0_2*eta_d0_2 + i_q0_2*eta_q0_2)/m_t2)


%--------------------------------------------------------------------------

% Condición de operación del sistema no lineal.
disp('La condición de operación del sistema no lineal es:');
fprintf('V_DC0 = %.4f\n', V_DC0);
fprintf('i_d0_1 = %.4f\n', i_d0_1);
fprintf('i_q0_1 = %.4f\n', i_q0_1);
fprintf('i_d0_2 = %.4f\n', i_d0_2);
fprintf('i_q0_2 = %.4f\n', i_q0_2);
fprintf('eta_d0_1 = %.4f\n', eta_d0_1);
fprintf('eta_q0_1 = %.4f\n', eta_q0_1);
fprintf('eta_d0_2 = %.4f\n', eta_d0_2);
fprintf('eta_q0_2 = %.4f\n', eta_q0_2);
fprintf('dV_DC_d = %.4f\n', dV_DC_d);

%--------------------------------------------------------------------------

% Controlador VSC1.
% Matrices principales del sistema lineal.
% Matriz A
A1 = [-R_1/L_1,  omega_1;
     -omega_1, -R_1/L_1];

% Matriz B
B1 = [-V_DC0/(2*L_1*m_t1),  0;
      0, -V_DC0/(2*L_1*m_t1)];

% Matriz C 
C1 = eye(2);

% Matriz D 
D1 = zeros(2);

% Acción integral.
% Matrices principales ampliadas del sistema linealizado.
Aa1 = [A1 zeros(2,2); -C1 zeros(2,2)]; 
Ba1 = [B1;-D1]; 
Ca1 = [C1 zeros(2,2)]; 
Da1 = D1; 

% Polos deseados del sistema en lazo cerrado.
Polos_deseados1 = [-2000 + 100i, -2000 - 100i, -1000 + 200i, -1000 - 200i];

% Verificacion de los polos deseados del sistema en lazo cerrado.
format shortg
disp('Polos deseados del sistema en lazo cerrado:');
disp(Polos_deseados1);
format

% Cálculo de la matriz de realimentación de estados Kca_place.
Kca_place1 = place(Aa1, Ba1, Polos_deseados1);

Kc1 = Kca_place1(:, 1:2); % Acción proporcional.
Ki1 = -Kca_place1(:, 3:4); % Acción integral.

format long e

% Verificacion del controlador.
disp('Matriz de ganancias Kc1:');
disp(Kc1);

% Verificacion del controlador.
disp('Matriz de ganancias Ki1:');
disp(Ki1);

% Polos del sistema en lazo cerrado.
Polos_modal1 = eig(Aa1-Ba1*Kca_place1);

% Verificacion de los polos del sistema en lazo cerrado.
format long e
disp('Polos del sistema en lazo cerrado:');
disp(Polos_modal1);
format


%--------------------------------------------------------------------------

% Controlador VSC2.
% Matrices principales del sistema lineal.
% Matriz A
A2 = [-R_2/L_2,  omega_2;
     -omega_2, -R_2/L_2];

% Matriz B
B2 = [-V_DC0/(2*L_2*m_t2),  0;
      0, -V_DC0/(2*L_2*m_t2)];

% Matriz C 
C2 = eye(2);

% Matriz D 
D2 = zeros(2);

% Acción integral.

% Matrices principales ampliadas del sistema linealizado.
Aa2 = [A2 zeros(2,2); -C2 zeros(2,2)]; 
Ba2 = [B2;-D2]; 
Ca2 = [C2 zeros(2,2)]; 
Da2 = D2; 

% Polos deseados del sistema en lazo cerrado.
Polos_deseados2 = [-2000 + 100i, -2000 - 100i, -1000 + 200i, -1000 - 200i];

% Verificacion de los polos deseados del sistema en lazo cerrado.
format shortg
disp('Polos deseados del sistema en lazo cerrado:');
disp(Polos_deseados2);
format

% Cálculo de la matriz de realimentación de estados Kca_place.
Kca_place2 = place(Aa2, Ba2, Polos_deseados2);

Kc2 = Kca_place2(:, 1:2); % Acción proporcional.
Ki2 = -Kca_place2(:, 3:4); % Acción integral.

format long e

% Verificacion del controlador.
disp('Matriz de ganancias Kc2:');
disp(Kc2);

% Verificacion del controlador.
disp('Matriz de ganancias Ki2:');
disp(Ki2);

% Polos del sistema en lazo cerrado.
Polos_modal2 = eig(Aa2-Ba2*Kca_place2);

% Verificacion de los polos del sistema en lazo cerrado.
format long e
disp('Polos del sistema en lazo cerrado:');
disp(Polos_modal2);
format

%--------------------------------------------------------------------------

% Diseño del filtro al set-point.
te = 1e-3; % Tiempo estacionario.
tau = te; % Constante de tiempo del filtro en segundos.

% Acción de control.
control = eye(2); % 1 ó 0.

%--------------------------------------------------------------------------

% Ruido en la medida de los sensores.
% Ruido sensores corrientes.
D_R_i = 5e-3;                 % Desviación tipica ruido i.
V_R_i = D_R_i^2;              % Varianza ruido i.

% Ruido salida controlador.
D_R_V = 2e-3;                 % Desviación tipica ruido voltaje.
V_R_V = D_R_V^2;              % Varianza ruido voltaje.

%--------------------------------------------------------------------------

% Desviaciones paramétricas.
% R_1:
Porc_Desv_R_1 = 0.05;             % Porcentaje desviacion.
Inst_Desv_R_1 = 0.1;              % Instante desviacion.

% L_1:
Porc_Desv_L_1 = 0.05;             % Porcentaje desviacion.
Inst_Desv_L_1 = 0.1;              % Instante desviacion.

% R_1:
Porc_Desv_R_2 = 0.05;             % Porcentaje desviacion.
Inst_Desv_R_2 = 0.1;              % Instante desviacion.

% L_1:
Porc_Desv_L_2 = 0.05;             % Porcentaje desviacion.
Inst_Desv_L_2 = 0.1;              % Instante desviacion.

% C_DC:
Porc_Desv_C_DC = 0.05;             % Porcentaje desviacion.
Inst_Desv_C_DC = 0.1;              % Instante desviacion.

%--------------------------------------------------------------------------

% Perturbacion externa.
I_p_i_d_1 = 0.02;         % Instante tiempo perturbación externa.
p_i_d_1 = 1.5;            % Magnitud perturbacion i_d_1.

I_p_i_q_1 = 0.02;         % Instante tiempo perturbación externa.
p_i_q_1 = -1.3;           % Magnitud perturbacion i_q_1.

I_p_i_d_2 = 0.02;         % Instante tiempo perturbación externa.
p_i_d_2 = 1.7;            % Magnitud perturbacion i_d_2.

I_p_i_q_2 = 0.02;         % Instante tiempo perturbación externa.
p_i_q_2 = -1.4;           % Magnitud perturbacion i_q_2.


%--------------------------------------------------------------------------
