Оптимизация PID для системы MIMO в MATLAB

Question

Я пишу код для оптимизации усиления PID системы MIMO с использованием метода наименьших квадратов в MATLAB. Завод, для которого я разрабатываю ПИД-регулятор, имеет 3 входа и 3 выхода, и для каждого входа я разрабатываю ПИД-контроллер. Вся система выглядит так:

%                   _____________________________________________________________________
%                  |                                                                     | 
%                  |                                                                     |  
%                  |                                                                     |
%                  |                      __________________________________             |  
%                  |    e              u |                                  |   y        |  
% r_phi,pitch,psi--|-->O--> PID array -->|              ssModel             |----------.-|------
%                  |   ^                 |                                  |          | |  
%                  |   | y               |__________________________________|          | |  
%                  |   |                                system                         | |  
%                  |   |                                                               | |   
%                  |   .---------------------------------------------------------------. |                                                      
%                  |_____________________________________________________________________|
%                                                        clSys
%    

Моя проблема в том, что логика оптимизации всегда возвращает оптимизированный коэффициент усиления, очень близкий к начальным значениям, по-видимому, он всегда находится в локальном минимуме.

Как я могу получить 9 выигрышей в глобальном масштабе?

Вот файл модели пространства состояний:

function [ssModel, Ts] = getSSmodel()
%% Define the parameter of our hardware
Jxx = 0.01; % appr. the inertia of a solid cube with the quad's mass (the best)
Jyy = 0.01;
Jzz = Jxx*2;

l   = 0.225;  % meters, the distance from the center of the propeller to the center of
              % the drone

Ts  = 0.0083; % the sampling time

b1  = l/Jxx;
b2  = l/Jyy;
b3  = 1/Jzz;

%% Define the continuous and discrete SS model
A = [0 1 0 0 0 0;
     0 0 0 0 0 0;
     0 0 0 1 0 0;
     0 0 0 0 0 0;
     0 0 0 0 0 1;
     0 0 0 0 0 0];

 B = [0  0  0;
      b1 0  0;
      0  0  0;
      0  b2 0;
      0  0  0;
      0  0  b3];

  C = [1 0 0 0 0 0;
       0 0 1 0 0 0;
       0 0 0 0 1 0];

  D = zeros(size(C,1),size(B,2));

  continuous_sys = ss(A,B,C,D);

  discrete_sys   = c2d(continuous_sys,Ts);

  ssModel     = discrete_sys;
end

Вот код, где я указываю функцию минимизации и функцию стоимости:

function costFunction = closedLoopPoles(pidGains, ssModel)

if (nargin < 2) 
    [ssModel, Ts]   =  getSSmodel();
end

%Define the PID controllers array
pidTF_array  = ones(3,3) * tf(0,1,Ts);

for i = 1:size(pidGains,1)
    pidTF_array(i,i)  = pid( pidGains(i,1), pidGains(i,2), pidGains(i,3),'IFormula','BackwardEuler','DFormula','BackwardEuler', 'Ts', Ts); 
end

% Assign the signals names and the systems structure
pidTF_array.u = 'e';          pidTF_array.y = 'u'; 
ssModel.u     = 'u';          ssModel.y     = 'y';
sum1          = sumblk('e = r - y',3);

system      = ssModel;

% Find the closed loop system
clSys               = connect(system,pidTF_array,sum1,'r','y');
clPoles            = pole( clSys );    
maxAbs_clPoles     = max(abs(clPoles));
sumAbs_clPoles     = sum(abs(clPoles));
sumnorm_clPoles    = sum(norm(clPoles));

costFunction       = sumnorm_clPoles;

end

И, наконец, вот функция оптимизации:

function pidGains_opt = optimizePID_minPole()

clear all
clc

pidGains_ini       = [180 4 20   160 4 20    150 4 10];
func               = @closedLoopPoles;
[pidGains_opt, ~]  = lsqnonlin(func,pidGains_ini)

end

Когда я запускаю код оптимизации, я получаю это в рабочей области:

 pidGains_ini =

   180     4    20
   160     4    20
   150     4    10

Warning: Trust-region-reflective algorithm requires at least as many equations as variables; using
Levenberg-Marquardt algorithm instead. 
> In lsqncommon at 77
  In lsqnonlin at 235
  In optimizePID_minPole at 11 

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the default value of the function tolerance.

<stopping criteria details>


pidGains_opt =

  179.8817    4.0022    3.1099
  159.8839    4.0025    3.0151
  149.9482    4.0022    1.5592