Using genetic algorithm (ga function) to generate a vector K considering a constraint

I'm trying to generate a vector K considering a constraint. My objective function is:

function k_bar = k_barra2(K)

A = [-0.5,0.4;-0.4,-0.5];
B = [0.0;2.5];
C = [0.0,2.6];
D = [0.0];
p = 5;

A = A - B * K;
C = C - D * K;

r = 1;
Id = eye(size(A,1));
Id_A = Id-A;
yss = (C*inv(Id_A)*B+D)*r;

[Y, X]=dstep(A,B,C,D,1,300);

Y = [Y; 0; 0];
eigen = eig(A);
for k=1:length(A)
    lambda = eigen(k);
    if((real(lambda) >= 0 && imag(lambda) == 0) && abs(lambda)<=1)
        isEigPos = 1;
    else
        isEigPos = 0;
        break;
    end
end
if(isEigPos == 1)
    if(yss > 0)
        infi = (yss - (yss * (p/100)));
        sup = (yss * (p/100) + yss);
    else
        sup = (yss - (yss * (p/100)));
        infi = (yss * (p/100) + yss);
    end

      i=1;
      while(~(Y(i)<=sup && Y(i)>=infi))
          i=i+1;
      end

      k_bar = i;
  else
      k=1;
      if(isEigPos == 1)
          Mp = yss;
          kp = k;
      else
          pre = Y(k);
          cur = Y(k+1);
          pos = Y(k+2);
          Mp = pre;
          kp = k;
          peak = pre;

          while((abs(Mp) <= abs(peak)) && (Mp ~= cur) && cur ~= Y(length(Y)))
              if((abs(cur) >= abs(pos)) && (abs(cur) >= abs(pre)))
                  peak = cur;
              end
              if ((Mp ~= peak)&&((yss>=0 && Y(k+1)>=0)||(yss<0 && Y(k+1)<0))&&(abs(peak) > abs(Mp)))
                  Mp = peak;
                  kp = k+1;
              end
              k = k+1;
              if k+2<=length(Y)
                  pre = cur;
                  cur = pos;
                  pos = Y(k+2);
              else
                  pre = Y(length(Y));
                  cur = Y(length(Y));
                  pos = Y(length(Y));
              end
          end
      end
      if Mp >= yss
          mp=Mp-yss;
      else
          mp=yss-Mp;
      end

      c_bar = (Mp-yss)/((max(abs(eig(A))))^(kp));
      x = abs((p*yss)/(100*c_bar));
      kss = log10(x)/log10(max(abs(eig(A))));
      k_bar = abs(ceil(kss))+size(A,1)

end
end

I want that the return of k_barra2() less or equal to a constant value ksr (could be, for example, ksr=5).

I tried to use the ga function as follows:

FitFcn = @k_barra2;
nvars = 2;
ksr = 5;
[K, fval] = ga(FitFcn, nvars, [1 1], ksr);

I know I'm wrong. It does not generate the right vector K. Anyone have an idea to solve that? Thanks in advance!

[Walter Roberson]

Do you mean that you want to stop optimization when the function result becomes as low as the threshold? As in you do not need the smallest practical value, just one that is Good Enough?

[Thiago Cavalcante]

Yes, It can be a K that is good enough, so satisfying k_barra2(K) <= ksr, where ksr can be 5. Thanks for the reply.

Answer 1

[Walter Roberson]

Use an output function that changes the state structure (first output, second input.) The documentation indicates,

    state — Structure containing information about the current generation. The State Structure describes the fields of state. To stop the iterations, set state.StopFlag to a nonempty character vector, such as 'y'.

[Thiago Cavalcante]

I did not get your point. Can you explain more, looking at my code above?

[Walter Roberson]

opt = optimoptions('ga', 'OutputFcn', @(options,state,flag) StopIfGoodEnough(options,state,flag,ksr));
A = []; b = [];
Aeq = []; beq = [];
lb = []; ub = [];
nonlcon = []; intcon = [];
[K, fval] = ga(FitFcn, nvars, A, b, Aeq, beq, lb, ub, nonlcon, intcon, opt);

Where

function [state, options, optchanged] = StopIfGoodEnough(options, state, flag, threshold)
  optchanged = false;
  if min(state.Best) <= threshold
    state.StopFlag = 'Good Enough';
  end

[Thiago Cavalcante]

Thanks for your reply. How can I consider that my K must satisfy to abs(eig(A-B*K)) == 1? Where A and B is inside k_barra2().

[Walter Roberson]

You should program that as a nonlinear constraint through the nonlcon parameter.

You should expect that ga will probably have a difficult time finding solutions to that.

[Thiago Cavalcante]

I tried to do what you suggested but I`m doing something wrong.

function [c,ceq]=stableSys(K)
A = [-0.5,0.4;-0.4,-0.5];
B = [0.0;2.5];
ceq = [];
  if abs(eig((A-B*K)))==1
    c=true
  else
    c=false
  end
end

And I replaced like:

[K, fval] = ga(FitFcn, nvars, A, b, Aeq, beq, lb, ub, @stableSys, intcon, opt);

[Walter Roberson]

Nonlinear equality constraints are considered to be satisfied if the value returned in ceq is sufficiently close to 0. It uses the value returned to hint on the direction to search. You are returning [] for your nonlinear equality constraint and are instead returning a value for the nonlinear inequality constraints.

function [c,ceq]=stableSys(K)
A = [-0.5,0.4;-0.4,-0.5];
B = [0.0;2.5];
c = [];
ceq = abs(eig((A-B*K)))-1;

[Thiago Cavalcante]

thanks for the reply Walter. The K returned by ga() did not work. It did not satisfy the 2 constraints (kbarra2()<=ksr and abs(eig((A-B*K)))==1*K)))==1).