Returning an additional value that is not part of the fitness function or objective for all population evaluation in GA

Can anyone please help by showing the code how can I return an additional value (Sigma) that is not part of my fitness function (Mean). I know how to return all the populations, scores but I'm not sure how to return this additional value (sigma). I see many have asked this question but no one showed in a code rather than just directing us to use the nested functions. Here is my code:

% the function to be optimized

 [objective]   = Optimization_Function(x,Pi,Pa,LogG,d,lifetime,Demand,BasePrice,HighPrice,...
  LowPrice,dis_rate_lamda,Geo,Wells_cost,Wells_rate,DStage,Operating_Fields,row,Amt,FOPEX,VOPEX)
.
.
.
.


expectedNPV  = mean(NPV_C);
sigma        = std(NPV_C); %%%%% I need this value for every population evaluation 
NewObjective = (expectedNPV)
objective    = - NewObjective; 

% my optimization code 
clear gaoutfunction
options          = optimoptions('ga','OutputFcn',@gaoutfunction,'UseParallel',true);
startTime        = tic;

fun              = @(x)Optimization_Function(x,Pi,Pa,LogNormal_G,d_cline,lifetime,Demand,BasePrice,...
    HighPrice,LowPrice,dis_rate_lamda,Geo,Wells_cost,Wells_rate,DStage,Operating_Fields,row,Amt,FOPEX,VOPEX);  

[xGA,fval]       = ga(fun,nvars,[],[],[],[],lowbond,upbond,[],[],options);
time_ga_parallel = toc(startTime);

record                = gaoutfunction();
gapopulationhistory   = vertcat(record.Population);
gabesthistory         = vertcat(record.Best);
gascorehistory        = vertcat(record.Score);
Results               = [gapopulationhistory gascorehistory];

% my output funciton which includes populations, scores but does not include sigma :(

function [state,options,optchanged] = gaoutfunction(options,state,flag)
    persistent state_record 
    if isempty(state_record)
      state_record = struct('Population', {}, 'Best', {}, 'Score', {});
    end
    if nargin == 0
      state = state_record;
      options = [];
      optchanged = [];
    else
      state_record(end+1) = struct('Population', state.Population, 'Best', state.Best', 'Score', state.Score);
      optchanged = false;
    end
end 

Can anyone please show me how to return sigma as I’m returning all the populations and their scores. Please note the score here is only the objective which is the mean. The sigma value is not part of the optimization but I need it to save the time rather than running the model again to evaluate it. Please help.

Answer 1

[Matt J]

function run_optimization

    allSigmas=[];

    [xGA,fval]   = ga(fun,nvars,[],[],[],[],lowbond,upbond,[],[],options);
    
    
     % the function to be optimized
    function  [objective]   = Optimization_Function(x,Pi,Pa,...)        .


                     ....
                         
           objective    = - NewObjective; 

          allSigmas(end+1)      = std(NPV_C); %%%%% I need this value for every population evaluation 
        
    end


end

[Yaser Khojah]

Dear Matt,

Thanks for your help but this does not work. First, how do you pass allSigmas to the main optimization_function? second how do you return allSigmas? somethings is missing.

[Matt J]

Why do you think something is missing? Did you try it?

There is no need to pass allSigmas back and forth beween Optimization_Function() and run_optimization(). Because Optimization_Function() is nested in run_optimization(), it shares access to allSigmas. Any change to allSigmas that you make in the workspace of Optimization_Function() is felt in the workspace of run_optimization() as well.

https://www.mathworks.com/help/matlab/matlab_prog/nested-functions.html#f4-73993

[Yaser Khojah]

Thanks Mat, the code did not did not return anything by implementing by your method. It was an empy cell. Would you please show what is missing as below:

%The optimizer
% %%%%%%%%%%%%%%%%%%%%%%%%%
allSigmas=[];
%%%%%%%%%%%%%%%%%%%%%%%%%%

fun              = @(x)Optimization_Function(x,Pi,Pa,LogNormal_G,d_cline,lifetime,Demand,BasePrice,...
    HighPrice,LowPrice,dis_rate_lamda,Geo,Wells_cost,Wells_rate,DStage,Operating_Fields,row,Amt,FOPEX,VOPEX);  

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

% the function to be optimized
 [objective]   = Optimization_Function(x,Pi,Pa,LogG,d,lifetime,Demand,BasePrice,HighPrice,...
  LowPrice,dis_rate_lamda,Geo,Wells_cost,Wells_rate,DStage,Operating_Fields,row,Amt,FOPEX,VOPEX)
.
.
.
.
expectedNPV  = mean(NPV_C);
sigma        = std(NPV_C); %%%%% I need this value for every population evaluation 
NewObjective = (expectedNPV)
objective    = - NewObjective; 

%%%%%%%%%%%%%%%%%%%%%%%%%%
allSigmas(end+1)  = sigma ;
%%%%%%%%%%%%%%%%%%%%%%%%%%

[Matt J]

I need to see the whole file in which

allSigmas=[]

resides.

[Yaser Khojah]

I put it in your format but I get only one value which is the initial value instead all of entire population (15200) as attached. Please check the first and last line where I reflected what you recommended as below. Thanks for your help!

function run_optimization

allSigmas=[];


%==========================================================================
%==========================================================================
%--------------------------------------------------------------------------
% *************NUMBER OF YEARS AND DISCOUNTE RATE WERE CAHANGE ************
%==========================================================================
%==========================================================================
%==========================================================================

% 
array_index        = '9'; 
% array_index        = getenv('PBS_ARRAY_INDEX');
iii                = str2num(array_index);




global drwaing 


Field       = 1;
SD_Range    = [0:4, 5:5:30];
% SD_Range    = 5;
array_index = size(SD_Range,2);
filename    = 'Fields_Input_Amt_OPT.json';

% for iii = 1:array_index

fprintf('%i\n', iii)


% the sampling was fixed here? please check
drwaing = 0; 


% filename = 'Fields_Input_Amt.json';  
[Pi,Pa,G,d,DStage,Wells_rate,Wells_cost,Geo,BasePrice,LowPrice,HighPrice,lifetime,Demand,dis_rate_lamda,DOFF,SOFF,On,Operating_Fields,Amt,FOPEX,VOPEX]= New_readinput_OPEX(filename);


% create the log normal for everyone
CE         = 0;
row        = 10; 
G_MC       = G/10^6; % convert G in place from MMSCF to TSCF
Mean       = G_MC; 
std_       = Mean .* CE; 
Variance   = std_.^2; 
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 

% create log-normal for 1 only
m            = G_MC(Field); 
sd           = SD_Range(iii);
v            = sd^2;
mu           = log((m^2)/sqrt(v+m^2));
sigma        = sqrt(log(v/(m^2)+1));
GIP(:,Field) = lognrnd(mu,sigma,[row,1]);
       
% put the unit back to MMSCFD
LogNormal_G      = GIP * 10^6;
rmax             = G .*d .* ones(size(LogNormal_G));
d_cline          = rmax ./ LogNormal_G; 

fnm = sprintf('GIP_%d.mat',SD_Range(iii));
save(fnm,'LogNormal_G');

fnm = sprintf('D_%d.mat',SD_Range(iii));
save(fnm,'d_cline');

% Update the pressure 
% pi_updated =  update_pressure(Pi, G, LogNormal_G, row);

% Run Optimizer

ub           = [3500, 2000, 1000, 3500, 2500,1000, 700, 2500];
lb           = [2000,  700,  400, 1000, 1000, 300, 300,  800];


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

clear gaoutfunction
options          = optimoptions('ga','OutputFcn',{@gaoutfunction},'UseParallel',true);
startTime        = tic;


fun              = @(x)Optimization_Function(x,Pi,Pa,LogNormal_G,d_cline,lifetime,Demand,BasePrice,...
    HighPrice,LowPrice,dis_rate_lamda,Geo,Wells_cost,Wells_rate,DStage,Operating_Fields,row,Amt,FOPEX,VOPEX);  

[xGA,fval]       = ga(fun,nvars,[],[],[],[],lowbond,upbond,[],[],options);
time_ga_parallel = toc(startTime);

record                = gaoutfunction();
gapopulationhistory   = vertcat(record.Population);
gabesthistory         = vertcat(record.Best);
gascorehistory        = vertcat(record.Score);
Results               = [gapopulationhistory gascorehistory];


fnm = sprintf('Optimization_Case_%d.mat',SD_Range(iii));
save(fnm);

fnm = sprintf('Results_%d.mat',SD_Range(iii));
save(fnm,'Results');


function [objective]   = Optimization_Function(x,Pi,Pa,LogG,d,lifetime,Demand,BasePrice,HighPrice,...
  LowPrice,dis_rate_lamda,Geo,Wells_cost,Wells_rate,DStage,Operating_Fields,row,Amt,FOPEX,VOPEX)




r_s          = x(1:Operating_Fields); 
t_0          = floor(x(Operating_Fields+1:end)); 
column       = size(r_s,2);
Supply       = zeros(lifetime+1,row);  
Revenue      = zeros(row,1);  
OPEX         = zeros(row,1);  
CapitalRisk  = zeros(row,1);
% F_off        = zeros(row,1);






[t_SS , t_A] = time_resevoir_Stochastic_Function(Pi,Pa,LogG,r_s,t_0,d,row,column);


% CAPEX = Total discoutned
% Dis_CAPEX = Discounted per field 
% CAPEX_Field = Captial per field without discounting

[CAPEX ,Dis_CAPEX, CAPEX_Field]  = OPT_CAPEX(r_s,t_0,Geo,Wells_cost,Wells_rate,dis_rate_lamda,DStage,row,column,Amt); 


for i = 1:row

[Supply_Cum, Field_No_Prod,Production_Profile_Cum,TimeSteps] = OPT_Supply(lifetime,t_0,t_SS(i,:),t_A(i,:),r_s,d);%^^^^^^
 Supply(:,i)            = Supply_Cum; 
[~,~,~,Reveune_Total,~] = New_Penalty_Average_Price_Function(Demand,Supply_Cum,BasePrice, HighPrice,LowPrice,dis_rate_lamda,lifetime);
Revenue(i)              = Reveune_Total; 


% The OPEX_ros is the total discounted rate 
[OPEX_row,~]   = OPT_OPEX(Production_Profile_Cum, CAPEX_Field(1,:),VOPEX, FOPEX,TimeSteps,lifetime,dis_rate_lamda);  %^^^^^^
OPEX(i)        = OPEX_row ; %^^^^^^
CAPEX_risk     = Dis_CAPEX(i,(Field_No_Prod));
CapitalRisk(i) = sum(CAPEX_risk,2);

%     if Field_No_Prod>0
%     F_off(i,1)     = Field_No_Prod; 
%     end 

end 

NPV_C        = (Revenue - CAPEX - OPEX); %^^^^^^
expectedNPV  = mean(NPV_C);
sigma        = std(NPV_C);
NewObjective = (expectedNPV + 0*sigma); 
objective    = - NewObjective;



allSigmas(end+1) = std(NPV_C); 



end 
end

allisgmas is only one value while my rest of populations are 15200 as attaced:

[Matt J]

Well, I can't run your code because not all required inputs have been provided, but below is a simplified example that should convince that you that the general technique does accumulate values into allSigma when properly implemented.

function test


   allSigmas=[];
   
   gaSimple(@Optimization_Function)
       
   allSigmas,
   

   
    function Optimization_Function(i)
        allSigmas(end+1)=i*10;
    end

end

function gaSimple(hfunc) 

   for i=1:10
       hfunc(i)
   end 

end

When run, the code should display,

allSigmas =

    10    20    30    40    50    60    70    80    90   100

[Yaser Khojah]

In this case yes you will return only a value which is only the initial value the optimizer tried. However, I will not be able to see all the stigmas that are related to every function evaluation. As you saw from the screenshot, I can see the whole 15000 populations and their scores. So, how can I see each sigma that is related to each function evaluation?. In your answer, the all sigma is only one value and I’m looking to find all stigmas are matching each population. One way to do that is to re run the main function by supplying all the population to get sigma but that takes a long time. I began to understand this is a limitation by MATLAB since no one has ever answered this question. You can try this in any problem to find if there is a way to return another value (not a single value; all the values that meet the population) that is not part of input nor the output of the optimizer. The thing I can print them but I cant return them.

[Matt J]

No, in the code I showed you a sigma value will be evaluated and stored every time OptimizationFunction is called. You can wrap another loop around gaSimple

   for j=1:3
   gaSimple(@Optimization_Function)
   end

to verify that this time you will get 30 values instead of 10,

   allSigmas =

  Columns 1 through 21

    10    20    30    40    50    60    70    80    90   100    10    20    30    40    50    60    70    80    90   100    10

  Columns 22 through 30

    20    30    40    50    60    70    80    90   100

[Yaser Khojah]

Dear Matt, thank you so much. It is finally working. I found the reason that stopped from work. My model was running in parallel and that makes it hard to store. So, switched the parallel off but if there is away please let me know. You just solved a year issue today. Thanks a lot. I gave up since no one was able to answer this question.

[Matt J]

I'm glad you solved it.

Regarding a workaround for UseParallel, the only solution I can think of right now is to switch on UseVectorized instead. Then, within your fitness function, implement your own explicit parfor loop over the members of the population:

function run_optimization


  allSigmas=[];

  opts = optimoptions('ga','UseVectorized','true');
  
  x=ga(@fitness,...,opts);

        function fitness(X)
        
        
          parfor i=1:size(X,1) 
        
              x=X(i,:); %current population member
              
              ...
                  
          
             y=  %operations on x
        
             sigma_vector(1,i)=std(y);
        
          end
        
          allSigmas(end+1,:)=sigma_vector;
          
        end

end

[Yaser Khojah]

Dear Matt, thanks for sharing this. I have tried and it is working. I never used parfor before. Do you think 'UseVectorized' is faster than using the parallel?

[Matt J]

UseVectorized can be faster than parfor-parallelization when the fitness of all population members can be computed all at once without a for-loop. I assumed that that was not the case for you.

[Yaser Khojah]

I agree, thanks :)