Aproximación de problemas de control. Cálculo de control óptimo. Resolución de la actividad propuesta.

Este artículo está en versión beta. El autor de este artículo no lo ha terminado todavía, por favor no lo edites hasta que elimine este mensaje.

// ARTICULO EN PREPARACIÓN

%
% Ejercicio sobre problemas inversos de control doctorado IMEIO
%

% Datos del problema

A = [ 0 1
     -1 1 ];
B = [ 1 0 ]';
t0 = 0;
T = 4;
Y0 = [ 1 1 ]';

% 1. Tenemos n == 2, m == 1

K = [ B A*B ];
rank ( K )

% rango ( K ) == m, cumple la condicion de Kalman y el sistema es controlable.

% 2. Solucion sin control por el método de Euler

N = 100;
h = ( T - t0 ) / N;
y = Y0;
y1 ( 1 ) = y ( 1 );
y2 ( 1 ) = y ( 2 );
for n = 1 : N
    y = y + h * A * y ;
    y1 ( n + 1 ) = y ( 1 );
    y2 ( n + 1 ) = y ( 2 );
end

% representación de la solución numerica:
figure ( 1 );
t = t0 : h : T;
plot ( t, y1, 'b')
hold on
plot ( t, y2, 'g')

% 3. Calcular un control u(t) para el cual y1(4) = 2; y2(4) = 4

YT = [ 2, 4 ]';

%  El control de minimo gasto energético viene dado por
%       u(t) := - B*(t) F*(t)^-1 F*(t) QT^-1 ( F(T)y0 - yT )
% con   F(t) = e^(tA) = 1+t^2A^2/2+t^3A^3/6+... matriz fundamental, F(t) = expm( t * A )
% y con QT = int_0_T ( F(t) B B* F*(t)* dt )

Ft = @( t ) expm ( t * A );
Ftt = @( t ) expm ( t * A' );
BB = B * B';
qt = @ ( t ) ( expm ( t * A ) * BB * expm ( t * A' ) );

% para la integral usamos la regla de Simpson compuesta
npanel = 32;
nn = 2 * npanel + 1;
hh = ( T - t0 ) / ( nn - 1 );
x = t0 : hh : T;
% definimos el gramiano a integrar para cada QT ( i, j )
for i = 1 : 1 : length ( x )
     q_t = qt ( x ( i ) );
     qt11 ( i ) = q_t ( 1, 1 );
     qt12 ( i ) = q_t ( 1, 2 );
     qt21 ( i ) = q_t ( 2, 1 ); 
     qt22 ( i ) = q_t ( 2, 2 );
end
%QT11 = (hh/3) * ( qt11( 1 ) + 4 * sum ( qt11( 2:2:nn-1 ) ) + 2 * sum ( qt11( 3:2:nn-2 ) ) + qt11( nn ) );
%QT12 = (hh/3) * ( qt12( 1 ) + 4 * sum ( qt12( 2:2:nn-1 ) ) + 2 * sum ( qt12( 3:2:nn-2 ) ) + qt12( nn ) );
%QT21 = (hh/3) * ( qt21( 1 ) + 4 * sum ( qt21( 2:2:nn-1 ) ) + 2 * sum ( qt21( 3:2:nn-2 ) ) + qt21( nn ) );
%QT22 = (hh/3) * ( qt22( 1 ) + 4 * sum ( qt22( 2:2:nn-1 ) ) + 2 * sum ( qt22( 3:2:nn-2 ) ) + qt22( nn ) );
%QT = [ QT11 QT12; QT21 QT22 ];

QT = 0;
for i = 1 : length ( t )
    Fti = Ft ( t ( i ) );
   QT = QT + Fti * BB * Fti' * h ;
end


%% integrate
u = zeros ( length ( t ), 1 );
%u = zeros ( length ( t ) );
FtT = Ft ( T );
for i = 1 : N
    Fti = Ft ( t ( i ) );
    u ( i ) = - B' * (((Fti)')^-1) * FtT' * QT^-1 * ( FtT * Y0 - YT );
end

figure ( 2 )
plot ( t, u, 'r' )

% 4. Solución  con  control por el método de Euler

N = 100;
h = ( T - t0 ) / N;
%yc = Y0 + B * u ( 1 );
yc = Y0 ;
yc1 ( 1 ) = yc ( 1 );
yc2 ( 1 ) = yc ( 2 );
for n = 1 : N
    yc = yc + h * ( A * yc + B * u ( n ) );
    yc1 ( n + 1 ) = yc ( 1 ); 
    yc2 ( n + 1 ) = yc ( 2 );
end

% representación de la solución numerica:

figure ( 1 )
plot ( t, yc1, '-- b' )
hold on
plot ( t, yc2, '-- g' )