% 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: Sistema de control no lineal.

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

% Subtítulo: MIMO.

% 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

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

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

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

% Subtítulo: IP Lazo P_2.

% Parámetros constantes para el diseño del controlador. 
t_e_P = 2e-2;          % Tiempo establecimiento (Segundos). 1e-2 -> 2e-2
tau_P = t_e_P/5;       % Constante de tiempo de un sistema de primer orden:
                       % t_e_P = 5*tau_P -> tau_P = t_e_P/5

% Constante de tiempo del sistema de control para la potencia activa:
% T_P (s) = 1/(((2+3*v_d_2*Kp_P)/(3*v_d_2*Ki_P)) s + 1)
% tau_P = ((2+3*v_d_2*Kp_P)/(3*v_d_2*Ki_P))
% Ki_P = ((2+3*v_d_2*Kp_P)/(3*v_d_2*tau_P))
% Calculo de Kp_P y Ki_P.
Kp_P = 0.0002; % 1 -> 0.0002
Ki_P = ((2+3*v_d_2*Kp_P)/(3*v_d_2*tau_P));

% Verificacion del controlador.
disp('Accion proporcional (Kp) del controlador de la potencia activa:');
disp(Kp_P);

% Verificacion del controlador.
disp('Accion integral (Ki) del controlador de la potencia activa:');
disp(Ki_P);

% Limites para la señal de control en amperios.
lim_i_d_2_upper = 75.07117;
lim_i_d_2_lower = -75.07117;

% Set-point P_2 en tres instantes diferentes.
P_2_1 = -3000;
P_2_2 = 3000;
P_2_3 = -8000;


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

% Subtítulo: IP Lazo Q_12.

% Parámetros constantes para el diseño del controlador.
t_e_Q = 2e-2;          % Tiempo establecimiento (Segundos). 1e-2 -> 2e-2 
tau_Q = t_e_Q/5;       % Constante de tiempo de un sistema de primer orden:
                       % t_e_Q = 5*tau_Q -> tau_Q = t_e_Q/5

% Constante de tiempo del sistema de control para la potencia reactiva:
% T_Q (s) = 1/(((3*v_d_12*Kp_Q-2)/(3*v_d_12*Ki_Q)) s + 1)
% tau_Q = ((3*v_d_12*Kp_Q-2)/(3*v_d_12*Ki_Q))
% Ki_Q = ((3*v_d_12*Kp_Q-2)/(3*v_d_12*tau_Q))
% Calculo de Kp_Q y Ki_Q.
Kp_Q = 0.0002; % 1 -> 0.0002
Ki_Q = ((3*v_d_1*Kp_Q-2)/(3*v_d_1*tau_Q));  % v_d_1 = v_d_2

% Verificacion del controlador.
disp('Accion proporcional (Kp) del controlador de la potencia reactiva:');
disp(Kp_Q);

% Verificacion del controlador.
disp('Accion integral (Ki) del controlador de la potencia reactiva:');
disp(Ki_Q);

% Limites para las señales de control en amperios.
% VSC1.
lim_i_q_1_upper = 393.72177;
lim_i_q_1_lower = -576.70686;

% VSC2.
lim_i_q_2_upper = 304.57446;
lim_i_q_2_lower = -446.13252;

% Set-point de Q_1 y Q_2 en tres instantes diferentes.
% VSC1.
% Q_1_1 = 0;
% Q_1_2 = -2500;
% Q_1_3 = -1500;
Q_1_1 = 0;
Q_1_2 = 0;
Q_1_3 = 0;

% VSC2.
% Q_2_1 = 0;
% Q_2_2 = -3000;
% Q_2_3 = -1000;
Q_2_1 = 0;
Q_2_2 = 0;
Q_2_3 = 0;


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

% Subtítulo: PI Lazo V_DC.

% Parámetros constantes para el diseño del controlador.
Xi = 1;                     % Factor amortiguamiento.
t_e_VDC = 3e-2;             % Tiempo establecimiento (Segundos). 1e-2 -> 3e-2
omega_n = 5.86/t_e_VDC;     % Frecuencia natural (Radianes/Segundos).

% Ecuación característica estándar de un sistema de segundo orden:
% s^2 + 2*Xi*omega_n s + omega_n^2 = 0
q_11 = 2*Xi*omega_n;
q_12 = omega_n^2;

% Ecuación característica del sistema en lazo cerrado:
% s^2 + (3*v_d_1)/(2*V_DC0*C_DC) * Kp s + (3*v_d_1)/(2*V_DC0*C_DC) * Ki = 0
q_21 = (3*v_d_1)/(2*V_DC0*C_DC);
q_22 = (3*v_d_1)/(2*V_DC0*C_DC);

% Calculo de Kp_DC y Ki_DC.
Kp_DC = q_11/q_21;
Ki_DC = q_12/q_22;

% Verificacion del controlador.
disp('Accion proporcional (Kp) del controlador de la tensión V_DC:');
disp(Kp_DC);

% Verificacion del controlador.
disp('Accion integral (Ki) del controlador de la tensión V_DC:');
disp(Ki_DC);

% Limites para la señal de control en amperios.
lim_i_d_1_upper = 75.07117;
lim_i_d_1_lower = -75.07117;

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

% Ruido en la medida de los sensores.
% Ruido sensor V_DC.
D_R_V = 30e-3;                 % Desviación tipica ruido voltaje.
V_R_V = D_R_V^2;               % Varianza ruido voltaje.

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

% Ruido salida controlador. Afecta demasiado en la potencia reactiva.
D_R_C = 2e-6;                 % Desviación tipica ruido voltaje. 2e-3;
V_R_C = 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.2;         % Instante tiempo perturbación externa.
p_i_d_1 = 2.5;            % Magnitud perturbacion i_d_1.

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

I_p_i_d_2 = 0.2;         % Instante tiempo perturbación externa.
p_i_d_2 = 2.7;            % Magnitud perturbacion i_d_2.

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


