My gamultiobj is not working but it works for the ga?

I have a function that its single objective was optimized by using ga optimizer in MATLAB. Everything is working. Now, I would like to solve the same function for two objectives by using multi objective ga optimizer (gamultiobj). When I run it, I get this massage “Reference to non-existent field 'Best'.”. Not sure what is wrong. Can anyone help please.

%% The funcitont is wokring when I use ga==== This is just to show the difference with single and multi objective
 function [objective]   = Stochastic_Model_Function_TOP(x,Pi,Pa,LogG,d,lifetime,Demand,BasePrice,HighPrice,...
  LowPrice,dis_rate_lamda,Geo,Wells_cost,Wells_rate,DStage,Operating_Fields,row)

options          = optimoptions('ga','OutputFcn',@gaoutfunction,'UseParallel',true);
startTime        = tic;
fun              = @(x)Stochastic_Model_Function_TOP(x,Pi,Pa,LogNormal_G,d,lifetime,Demand,BasePrice,...
    HighPrice,LowPrice,dis_rate_lamda,Geo,Wells_cost,Wells_rate,DStage,Operating_Fields,row);  
time_ga_parallel = toc(startTime)
record                = gaoutfunction();
gapopulationhistory   = vertcat(record.Population);
gabesthistory         = vertcat(record.Best);
gascorehistory        = vertcat(record.Score);

The code below is not working and not sure why?

%% Now I would like to use on Stochastic_Model_Function_TOP so I fixed the function as the following 

% inside the main funciton i added the objetives as below
[objective]   = Stochastic_Model_Function_TOP(x,Pi,Pa,LogG,d,lifetime,Demand,BasePrice,HighPrice,...
  LowPrice,dis_rate_lamda,Geo,Wells_cost,Wells_rate,DStage,Operating_Fields,row)

objective(1) = - expectedNPV;
objective(2) = sigma;

% now to optimize the funciton i added the below but it is not working. 

options    = optimoptions('gamultiobj','OutputFcn',@gaoutfunction,'UseParallel',true);
startTime  = tic;     
fun        = @(x)Stochastic_Model_Function_TOP(x,Pi,Pa,LogNormal_G,d,lifetime,Demand,BasePrice,HighPrice,...
    LowPrice,dis_rate_lamda,Geo,Wells_cost,Wells_rate,DStage,Operating_Fields,row);

time_ga_parallel = toc(startTime)
record                = gaoutfunction();
gapopulationhistory   = vertcat(record.Population);
gabesthistory         = vertcat(record.Best);
gascorehistory        = vertcat(record.Score);
         

      

[Alan Weiss]

Please give the ENTIRE error message, everything in red. And give your gamultiobj call.

Alan Weiss

MATLAB mathematical toolbox documentation

[Yaser Khojah]

Dear Alan thank you for your response and sorry for my late since I left the offfice.

This is the error massage I get:

Reference to non-existent field 'Best'.

Error in gaoutfunction (line 13)

state_record(end+1) = struct('Population', state.Population, 'Best', state.Best', 'Score',

state.Score);

Error in gaoutput (line 39)

[state,optnew,changed] = feval(functions{i},options.OutputPlotFcnOptions,state,flag,args{i}{:});

Error in gamultiobjsolve (line 13)

[state,options] = gaoutput(FitnessFcn,options,state,currentState);

Error in gamultiobj (line 303)

[x,fval,exitFlag,output,population,scores] = gamultiobjsolve(FitnessFcn,nvars, ...

Error in MeanVariance_OPT (line 71)

[xGA,fval] = gamultiobj(fun,nvars,[],[],[],[],lowbond,upbond,options);

and here is my code:

global drwaing 

drwaing = 0; 

filename = 'Fields_Input.json';  
[Pi,Pa,G,d,DStage,Wells_rate,Wells_cost,Geo,BasePrice,LowPrice,HighPrice,lifetime,Demand,dis_rate_lamda,DOFF,SOFF,On,Operating_Fields]= New_readinput(filename);


Persent_v = .25;
xGA_v     = zeros(length(Persent_v), Operating_Fields*2);
fval_v    = zeros(1,length(Persent_v));


for i = 1:length(Persent_v)

        Persent    = Persent_v(i); 
        row        = 1000; 
        G_MC       = G/10^6; % convert G in place from MMSCF to TSCF
        Mean       = G_MC; 
        Variance   = G_MC .* Persent;      
        
        mu      = log((Mean.^2)./sqrt(Variance+Mean.^2));
        sigma   = sqrt(log(Variance./(Mean.^2)+1));
        GIP     = zeros(row,length(G_MC));
        

        for index = 1:length(G_MC)
            
            R   = lognrnd(mu(index),sigma(index),[row,1]);
         
            
            GIP(:,index) = R; 

        end 
        
        LogNormal_G     = GIP * 10^6; % put the unit back to MMSCFD

        % Run Optimizer 
        ub           = [3835.61643835617,2301.36986301370,887.671232876712,4109.58904109589,2958.90410958904,1109.58904109589,863.013698630137,2739.72602739726];
        lb           = [1095.89041095890,605.623648161500,317.025440313112,821.917808219178,845.401174168298,277.397260273973,246.575342465753,608.828006088280];

        nvars         = Operating_Fields*2;
        MaxProduction = ub *365; 
        MinProduction = lb *365;
        upbond        = [MaxProduction  lifetime* ones(1,Operating_Fields)]; 
        lowbond       = [MinProduction  zeros(1,Operating_Fields)]; 
        
        clear gaoutfunction

        options          = optimoptions('gamultiobj','OutputFcn',@gaoutfunction,'UseParallel',true);

        startTime        = tic;
       
        fun              = @(x)Stochastic_Model_Function_TOP(x,Pi,Pa,LogNormal_G,d,lifetime,Demand,BasePrice,HighPrice,LowPrice,dis_rate_lamda,Geo,Wells_cost,Wells_rate,DStage,Operating_Fields,row);  
        
       [xGA,fval]        = gamultiobj(fun,nvars,[],[],[],[],lowbond,upbond,options);
       time_ga_parallel  = toc(startTime);
       
       record                = gaoutfunction();
       gapopulationhistory   = vertcat(record.Population);
       gabesthistory         = vertcat(record.Best);
       gascorehistory        = vertcat(record.Score);

end