Why does this code give me an error?

I want to run the Wild Horse Optimizer (WHO) with my fitness function. It will give me an estimated solution "gBest". Then I want to supply that estimated solution to the local builtin optimizer "sqp" while providing an initial guess as a starting point. My code is below:

clear;clc
u=[-33 33];% desired vector
dim=length(u);
lb=-90*ones(1,dim);
ub= 90*ones(1,dim);

Noise=0;

% Define the objective function
objectiveFunctionForMetaheuristic = @(x) MIMOfunNoise(x,u,Noise) + penaltyTerm(x, u);
objectiveFunctionForLocal = @(x) MIMOfunNoise(x, u, Noise) + penaltyTerm(x, u);

% Initial guess 
initialGuess = u;

% Set options for fmincon (SQP solver)
options = optimoptions('fmincon', 'Algorithm', 'sqp', 'Display', 'off');

Runs=30;

% Pre-allocation for WHO algorithm
one=zeros(Runs,1);
time=zeros(Runs,1);
temp=zeros(Runs,dim);
two=zeros(Runs,dim);

% Pre-allocation for SQP algorithm
oneSQP=zeros(Runs,1);
timeSQP=zeros(Runs,1);
tempSQP=zeros(Runs,dim);
twoSQP=zeros(Runs,dim);

nn=0;
for n=1:Runs %------------(2)
    nn=nn+1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Call WHO 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[time,gBest,gBestScore]=WHO(N,Max_iter,lb,ub,dim,objectiveFunction);
[time,gBest,gBestScore]=WHO(30,500,lb,ub,dim,objectiveFunctionForMetaheuristic);
one(nn)=gBestScore;
temp(nn,:)=gBest;
time1(nn)=time;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Swapping
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[~, ix] = sort(u); 
[~, ix1(ix)] = sort(temp(nn,:));  
two(nn,:) = temp(nn,ix1); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Call the SQP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic;
%[bestX, fmin] = fmincon(objectiveFunctionForLocal, gBest, [], [], [], [], lb, ub, [], options);
objectiveFunctionForLocal(gBest)
[c,ceq] = nonlinear_constraint(gBest)
[bestX, fmin] = fmincon(objectiveFunctionForLocal, gBest, [], [], [], [], lb, ub, @(x) nonlinear_constraint(x), options, initialGuess);
timeSQP(nn)=toc;
oneSQP(nn)=fmin;
tempSQP(nn,:)=bestX;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Swapping
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[~, ix] = sort(u); % u is my desired vector
[~, ix1(ix)] = sort(tempSQP(nn,:));  
twoSQP(nn,:) = tempSQP(nn,ix1); 
end

ans = 1.2321

Unrecognized function or variable 'nonlinear_constraint'.

%%%%%%%%%%%%%%%%%%%
% Choose best
%%%%%%%%%%%%%%%%%%%
[one1 ind]=sort(one,'descend');
[fitness,ind1]=min(one1);
two1=two(ind1,:);

[one1SQP indSQP]=sort(oneSQP','descend');
[fitnessSQP,ind1SQP]=min(one1SQP);
two1SQP=twoSQP(ind1SQP,:);


% Display results
format compact
fprintf('Global_Sol  : %s\n', num2str(round(two1, 4),'%.4f '));
fprintf('Local_Sol   : %s\n', num2str(round(two1SQP,4),'%.4f '));
fprintf('Desired     : %s\n', num2str(u,'%.4f '));
fprintf('Global_fmin : %f\n', fitness);
fprintf('Local_fmin  : %f\n', fitnessSQP);

%%%%%%%%%%%%%%%%%%
% Save workspace data
%%%%%%%%%%%%%%%%%
%save 2sn0dB

function e=MIMOfunNoise(b,u,Noise)
M=6; % const1
N=6; % const2
K=length(u);%const3
d = 0.5; % const4
vec = @(MAT) MAT(:); 
vecH = @(MAT) MAT(:).'; 
steerVecT = @(ang) exp(1j*2*pi*d*(0:M-1).'*sin(vecH(ang))); 
steerVecR = @(ang) exp(1j*2*pi*d*(0:N-1).'*sin(vecH(ang))); 
%%%%%%%%%%%%%%%%%%%%
% Swapping vector b
%%%%%%%%%%%%%%%%%%%%
[~, ix] = sort(u); 
[~, ix1(ix)] = sort(b);  
b = b(ix1);
A = ones(K, 1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculation of yo and ye
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yo = yMatTR(deg2rad(u), steerVecT, steerVecR); 
yo=awgn(yo,Noise);
ye = yMatTR(deg2rad(b), steerVecT, steerVecR); 
%%%%%%%%%%%%%%%%%%
% MSE
%%%%%%%%%%%%%%%%%%
e=norm(yo-ye).^2/(M*N);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% User defined function called within MIMOfunNoise
%%%%%%%%%%%%%%%%%%%%%%%%5%%%
function y = yMatTR(targetAngle, steerVecT, steerVecR)
   steerA = steerVecT(targetAngle); 
   steerB = steerVecR(targetAngle); 
   y=sum( steerA.*permute(steerB,[3,2,1]) ,2);
   y=y(:);   
end

function penalty = penaltyTerm(x, desiredResult)
    % Define a penalty term based on the deviation from the desired result
    penalty = 100 * sum((x - desiredResult).^2);
end

function [X]=initialization(N,dim,up,down)

if size(up,1)==1
    X=rand(N,dim).*(up-down)+down;
end
if size(up,1)>1
    for i=1:dim
        high=up(i);low=down(i);
        X(:,i)=rand(1,N).*(high-low)+low;
    end
end
end

function Stallion=exchange(Stallion)
nStallion=length(Stallion);

for i=1:nStallion 
   [value,index]=min([Stallion(i).group.cost]);
  
   if value<Stallion(i).cost
       
       bestgroup=Stallion(i).group(index);
       
       Stallion(i).group(index).pos=Stallion(i).pos;
       Stallion(i).group(index).cost=Stallion(i).cost;
       
       
       Stallion(i).pos=bestgroup.pos;
       Stallion(i).cost=bestgroup.cost;
   end   
end
end
% Developed in MATLAB R2017b
% Source codes demo version 1.0
% _____________________________________________________

%  Author, inventor and programmer: Iraj Naruei and Farshid Keynia,
%  e-Mail: [e-mail] , [e-mail]
% _____________________________________________________
%  Co-author and Advisor: Farshid Keynia
%
%         e-Mail: [e-mail]
% _____________________________________________________
%  Co-authors: Amir Sabbagh Molahoseyni
%        
%         e-Mail: [e-mail]
% _____________________________________________________
% You can find the Wild Horse Optimizer code at 
% _____________________________________________________
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Max_iter: maximum iterations, N: populatoin size, Convergence_curve: Convergence curve
%function [Convergence_curve,gBest,gBestScore]=WHO(N,Max_iter,lb,ub,dim,fobj)
function [time,gBest,gBestScore]=WHO(N,Max_iter,lb,ub,dim,fobj)
tic; % By Me
if size(ub,1)==1
    ub=ones(1,dim).*ub;
    lb=ones(1,dim).*lb;
end
 PS=0.2;     % Stallions Percentage 
 PC=0.13;    % Crossover Percentage
 NStallion=ceil(PS*N); % number Stallion
 Nfoal=N-NStallion;

Convergence_curve = zeros(1,Max_iter);
gBest=zeros(1,dim);
 gBestScore=inf;
 
 %create initial population
 empty.pos=[];
empty.cost=[];

group=repmat(empty,Nfoal,1);

for i=1:Nfoal 
   group(i).pos=lb+rand(1,dim).*(ub-lb);
   group(i).cost=fobj(group(i).pos);
end

Stallion=repmat(empty,NStallion,1);

for i=1:NStallion 
   Stallion(i).pos=lb+rand(1,dim).*(ub-lb);
   Stallion(i).cost=fobj(Stallion(i).pos);
 
end
  ngroup=length(group);
  a=randperm(ngroup);
  group=group(a);

i=0;
k=1;
for j=1:ngroup
i=i+1;    
  Stallion(i).group(k)=group(j);  
if i==NStallion
    i=0;
    k=k+1;
end
end
Stallion=exchange(Stallion);
[value,index]=min([Stallion.cost]);

WH=Stallion(index); % global
   gBest=WH.pos;
   gBestScore=WH.cost;

Convergence_curve(1)=WH.cost;
l=2; % Loop counter
while l<Max_iter+1
TDR=1-l*((1)/Max_iter);

for i=1:NStallion
    
   ngroup=length(Stallion(i).group);
    [~,index]=sort([Stallion(i).group.cost]);
    Stallion(i).group=Stallion(i).group(index);
   
   for j=1:ngroup
    
    if rand>PC
              z=rand(1,dim)<TDR;
            r1=rand;
            r2=rand(1,dim);
            idx=(z==0);
            r3=r1.*idx+r2.*~idx;
           rr=-2+4*r3;

           Stallion(i).group(j).pos= 2*r3.*cos(2*pi*rr).*(Stallion(i).pos-Stallion(i).group(j).pos)+(Stallion(i).pos);
    else
    A=randperm(NStallion);
    A(A==i)=[];
    a=A(1);
    c=A(2);
%     B=randperm(ngroup);
%     BB=randperm(ngroup);
%     b1=B(1);b2=BB(1);
    x1=Stallion(c).group(end).pos;
    x2=Stallion(a).group(end).pos;

       y1=(x1+x2)/2;   % Crossover

    Stallion(i).group(j).pos=y1;
    end
   

    Stallion(i).group(j).pos=min(Stallion(i).group(j).pos,ub);
    Stallion(i).group(j).pos=max(Stallion(i).group(j).pos,lb);
    
    Stallion(i).group(j).cost=fobj(Stallion(i).group(j).pos);
    

    
   end
    
% end
% 
% for i=1:NStallion

    R=rand;
%     z=rand(1,dim)<TDR;
%             r1=rand;
%             r2=rand(1,dim);
%             idx=(z==0);
%             r3=r1.*idx+r2.*~idx;
%             rr=-2+4*r3;

       if R<0.5
        k= 2*r3.*cos(2*pi*rr).*(WH.pos-(Stallion(i).pos))+WH.pos;
        else
        k= 2*r3.*cos(2*pi*rr).*(WH.pos-(Stallion(i).pos))-WH.pos;
       end

    k=min(k,ub);
    k=max(k,lb);
    fk=fobj(k);
    if fk<Stallion(i).cost
      Stallion(i).pos  =k;
      Stallion(i).cost=fk;
    end
end
    Stallion=exchange(Stallion);
     [value,index]=min([Stallion.cost]);
   if value<WH.cost
       WH=Stallion(index);
   end
   gBest=WH.pos;
   gBestScore=WH.cost;
    Convergence_curve(l)=WH.cost;
    l = l + 1;
    time=toc; % By Me
end
end

The WHO and its supported files are in the attachment. But when I run the above code, it gives me an error like this:

Error using abcd>@(x)MIMOfunNoise(x,u,Noise)+penaltyTerm(x,u)

Too many input arguments.

Error in fmincon (line 563)

initVals.f = feval(funfcn{3},X,varargin{:});

Error in abcd (line 61)

[bestX, fmin] = fmincon(objectiveFunctionForLocal, gBest, [], [], [], [], lb, ub, @(x) nonlinear_constraint(x), options, initialGuess);

Caused by:

Failure in initial objective function evaluation. FMINCON

cannot continue.

>>

[Torsten]

Your constraint function "nonlinear_constraint" is missing.

[Sadiq Akbar]

@Torsten

But what is meant by nonlinear_contraint function and hiw will we define it?

[Torsten]

You call "fmincon" and tell him that you want to supply a constraint function "nonlinear_constraint" that "fmincon" will search for. Where is this function ?

[bestX, fmin] = fmincon(objectiveFunctionForLocal, gBest, [], [], [], [], lb, ub, @(x) nonlinear_constraint(x), options, initialGuess);

And what is the "initialGuess" argument for if you want to supply gBest as initial guess ?

[Sam Chak]

Hi @Sadiq Akbar,

It appears that the code was originally authored by Iraj Naruei and Farshid Keynia. Could you please clarify which part of the code you have modified to suit your application? It is important to establish your understanding of the entire code and the extent of your modifications. Without this clarity, it may be challenging to effectively communicate and address programming issues.

Based on logical deduction, I suspect that you may have modified or supplied your own Cost Function, MIMOfunNoise(), but without including constraints in your optimization problem. In the authors' demo, the original fmincon() optimization problem likely required the nonlinear_constraint(x) function. If you run the code without addressing this incompatibility issue, an error is likely to occur.

It is unlikely that any designer who created the Cost Function would not understand whether or not there are constraints placed on the design. Could you please write out the optimization problem in mathematical form? By expressing the problem mathematically, you can later translate it into MATLAB code for implementation.

[Sadiq Akbar]

@Torsten

Thanks a lot fro your kind response. As I am not expert so much in Matlab, so I don't know what is nonlinear_contraint function? But the bounds of my variables are already defined to be :

lb=-90*ones(1,dim);
ub= 90*ones(1,dim);

Can you guide me what is meant by nonlinear constraint function. If you guide me, then I will try to define it but actually I don't know much about this function. Earlier when I used to use the GA tool, there also, I usually left it blank. But since there is no GA tool in my MATLAB now, so I cannot use the GA tool now. That's why I feel much difficulty now.

[Sadiq Akbar]

@Torsten

Actually I have downloaded the WHO algorithm from the Mathworks site. And the MIMOfunNoise() is my fitness function and I took the code of sqp() from the google but now I don't know the source where I took it from.

However as far as I know there is no nonlinear constraint function with my fitness function. If we remove it, how will we run it?

[Sam Chak]

@Sadiq Akbar, you should check and compare the fmincon() function in your code with the proper syntax of fmincon as shown below:

[bestX, fmin] = fmincon(objectiveFunctionForLocal, ...  % cost function
    gBest, ...                                          % initial point
    [], ...                                             % A
    [], ...                                             % b
    [], ...                                             % Aeq
    [], ...                                             % beq
    lb, ...                                             % lower bound
    ub, ...                                             % upper bound
    @(x) nonlinear_constraint(x), ...                   % nonlinear inequalities
    options, ...                                        % optimoptions()
    initialGuess);                                      % ???

[Rena Berman]

(Answers Dev) Restored edit

Answer 1

[Sam Chak]

Hi @Sadiq Akbar

If nonlinear inequality constraints are not available, you can disable them in the code. It appears that the updated code now returns some results. Please verify and check.

% [c,ceq] = nonlinear_constraint(gBest)
[bestX, fmin] = fmincon(objectiveFunctionForLocal, gBest, [], [], [], [], lb, ub, [], options)

bestX = 1×2

-32.9928 32.9909

fmin = 0.4869

bestX = 1×2

-32.9921 33.0007

fmin = 0.7671

bestX = 1×2

-32.9935 33.0065

fmin = 0.5644

bestX = 1×2

-32.9949 32.9937

fmin = 0.6559

bestX = 1×2

-33.0059 33.0040

fmin = 0.7550

bestX = 1×2

-32.9730 32.9957

fmin = 0.5571

bestX = 1×2

-32.9869 33.0073

fmin = 0.5580

bestX = 1×2

-33.0045 33.0003

fmin = 0.6424

bestX = 1×2

-32.9803 33.0095

fmin = 0.5441

bestX = 1×2

-33.0141 33.0049

fmin = 0.5860

bestX = 1×2

-33.0051 33.0009

fmin = 0.5768

bestX = 1×2

-33.0024 33.0039

fmin = 0.4507

bestX = 1×2

-33.0093 33.0032

fmin = 0.5387

bestX = 1×2

-32.9842 32.9942

fmin = 0.6914

bestX = 1×2

-32.9965 32.9991

fmin = 0.5265

bestX = 1×2

-32.9984 32.9997

fmin = 0.6837

bestX = 1×2

-32.9989 33.0143

fmin = 0.4914

bestX = 1×2

-32.9944 32.9951

fmin = 0.5838

bestX = 1×2

-33.0080 33.0096

fmin = 0.6222

bestX = 1×2

-32.9865 33.0004

fmin = 0.6634

bestX = 1×2

-33.0073 33.0058

fmin = 0.7632

bestX = 1×2

-33.0026 33.0117

fmin = 0.5873

bestX = 1×2

-32.9913 32.9825

fmin = 0.7897

bestX = 1×2

-33.0031 32.9977

fmin = 0.5744

bestX = 1×2

-33.0043 32.9873

fmin = 0.7019

bestX = 1×2

-33.0017 33.0097

fmin = 0.5797

bestX = 1×2

-33.0032 32.9998

fmin = 0.6127

bestX = 1×2

-33.0016 32.9977

fmin = 0.5354

bestX = 1×2

-33.0086 33.0055

fmin = 0.5288

bestX = 1×2

-32.9955 32.9920

fmin = 0.6432

%%%%%%%%%%%%%%%%%%%
% Choose best
%%%%%%%%%%%%%%%%%%%
[one1 ind]=sort(one,'descend');
[fitness,ind1]=min(one1);
two1=two(ind1,:);

[one1SQP indSQP]=sort(oneSQP','descend');
[fitnessSQP,ind1SQP]=min(one1SQP);
two1SQP=twoSQP(ind1SQP,:);


% Display results
format compact
fprintf('Global_Sol  : %s\n', num2str(round(two1, 4),'%.4f '));

Global_Sol  : -32.9952  32.9930

fprintf('Local_Sol   : %s\n', num2str(round(two1SQP,4),'%.4f '));

Local_Sol   : -32.9955  32.9920

fprintf('Desired     : %s\n', num2str(u,'%.4f '));

Desired     : -33.0000  33.0000

fprintf('Global_fmin : %f\n', fitness);

Global_fmin : 0.269830

fprintf('Local_fmin  : %f\n', fitnessSQP);

Local_fmin  : 0.450739

[Sadiq Akbar]

@Sam Chak

Thanks a lot for your guiadance. Yes there is no nonlinear constraint in my problem. But now again there is a problem. In earlier times, when I used the GA tool, and I used to run my code with GA, SQP and then GA_SQP. At that time the estimated solution by GA wa less accurate, likewise the SQP estimation was more accurate but the estimation by GA_SQP was the best. I mean when I gave the estimated solution of Ga to SQP as an input, then the SQP would fine tune that estimated solution further and therefore the estimated solution by the hybrid GA_SQP would be the best at that time but here it is not so. Why?

[Sam Chak]

You're welcome, @Sadiq Akbar. I simply fixed the fmincon syntax in the optimization code to get it running. As for why the GA-SQP hybrid is expected to perform better in your problem, I cannot provide a definitive answer at the moment.

Answer 2

[Walter Roberson]

[bestX, fmin] = fmincon(objectiveFunctionForLocal, gBest, [], [], [], [], lb, ub, @(x) nonlinear_constraint(x), options, initialGuess);

fmincon() does not expect an initialGuess parameter. It sees the initialGuess as an extra parameter. What it does with that extra parameter is no longer documented: it passes that extra parameter as an extra parameter to objectiveFunctionForLocal and as an extra parameter to @(x) nonlinear_constraint(x) too. But

objectiveFunctionForLocal = @(x) MIMOfunNoise(x, u, Noise) + penaltyTerm(x, u);

expects only a single parameter, and so gives an error when the extra parameter intialGuess is passed to it.

[Sadiq Akbar]

@Walter Roberson

Thanks a lot for your kind guiadance. Yes, I removed the nonlinear constraint function and it works but the problem is that when I run my fitness function MIMOfunNoise() with GA, it gives me a best estimated vector, then I give that best estimated vetor to the SQP so that SQP further fine tune it( as I used to do before) but it is not so now i.e., I mean the SQP gives even a worse solution and fitness than GA. I don't know why and where am I doing the mistake?