PID gains optimization for a MIMO system in MATLAB

Question

I am writing a code to optimize the PID gains of a MIMO system using least-squares in MATLAB. The plant I am designing the PID controller for has 3 inputs and 3 outputs, and for each input I am designing a PID controller. The whole system looks like this:

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

My problem is that the optimization logic always returns the optimized gains very close to the initial values, apparently it always fall in a local minimum.

How can I get the 9 gains globally minimized?

Here is the state space model file:

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

Here is the code where I specify the function to minimize and the cost function:

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

And finally here is the optimization function:

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

When I run the optimization code, I get this in the workspace:

 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