gamultiobj optimizer is not working and I get the following error massage "Reference to non-existent field 'Best'"

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). However, the model does not run and I ge the error massage as below:

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);
       

My script code as below:

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));

% the for loop since I will later have more than a persent_V value

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 

my main function to be optimized as below:

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)

%Cant include everything but I'm trying to show how my objective is linked as below

NPV_C        = (Revenue - CAPEX);
expectedNPV  = mean(NPV_C);
sigma        = std(NPV_C);

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

end