How to solve multi-objective optimization problem of a single variable that can take only integer values?

Hi,

I am trying to solve a multi-objective optimization problem where I only have a single variable that can take integer values. I tried the code provided by the Mathworks support team for imposing integer constraints which is written for two variables. I was able to use the functions provided in exampleIntGAMULTIOBJ.zip: @int_pop, @int_mutation, @int_crossoverarithmetic for solving problem with two variables. However, it doesn't work for the case of a single variable. What changes should I make to the functions provided to be able to make it work for single variable? Below is the link for the functions

https://www.mathworks.com/matlabcentral/answers/103369-is-it-possible-to-solve-a-mixed-integer-multi-objective-optimization-problem-using-global-optimizati

Thanks.

[Jeff Miller]

If there is only one integer variable to be adjusted, why not just try different values (e.g., with a for-loop across the acceptable range) and see which one works best?

[Vamshi Krishna Gudipati]

The range is huge and the computation cost for each loop is quite high.

[Jeff Miller]

Do you have any idea of the shape of the function that you are trying to optimize? E.g. what are its values for example integer parameter values, e.g., i=10,20,50,100,500,1000,10000,100000 (I'm just guessing about the range)? If the function is monotonic on both sides of the optimal point, then you should be able to narrow in on that point with something like a bisection method.

[Sachin Meena]

Did Jeff's suggestion work for you? If not, could you explain why not?

[Vamshi Krishna Gudipati]

Jeff, Thanks for your response. The function I am dealing with is non-linear and I don't have a clear idea of the shape of the function. Actually, I am more interested in trying to make the example code provided by Mathworks staff work for a single variable.

[Vamshi Krishna Gudipati]

Sachin,

I am interested to make the code provided in the above link to work for a single variable. I tried several possibilities but I keep on getting the following error. If you could help that would be great. Index exceeds matrix dimensions.

Error in objAndConVectorizer (line 32) if isFeas(i)

Error in gamultiobjMakeState (line 225) [nextScore,nextC,nextCeq,nextIsFeas] = objAndConVectorizer(nextPopulation, ...

Error in gamultiobjsolve (line 8) state = gamultiobjMakeState(GenomeLength,FitnessFcn,ConstrFcn,output.problemtype,options);

Error in gamultiobj (line 274) [x,fval,exitFlag,output,population,scores] = gamultiobjsolve(FitnessFcn,nvars, ...

Error in Pareto_Optimization_RL_IC_HB (line 30) [x, f, exitflag, output, population, score] = gamultiobj(FitnessFunction,NOV,A,b,Aeq,beq,lb,ub,options);

Caused by: Failure in user-supplied fitness function evaluation. Cannot continue.

Thanks, Vamshi.

[Sachin Meena]

Can you share your fitness function with sufficient comments?

[Cesar Antonio Lopez Segura]

Hi Vamshi,

Can you share the code of your objective function? I can test it with Matlab opens source ga code.

[Vamshi Krishna Gudipati]

function y = HB_RegionalLoss_InitialCost(x)

[sec_prop,~] = xlsread('AISC_Shapes_Database.xls','Database v14.1');
Upload1 = importdata('PGA_short.txt');  
Upload2 = importdata('Annual_exceedance_rate_short.txt');
PGA = Upload1.data;   % Discretized PGA values (same PGA range that is used to develop neural network )
ARE = Upload2.data;   % Annual rate of exceedance

% Probability of a pga value
PPGA(1) = 1-(1-exp(-0.528973));  % 0.528973 corresponds to 0.005 pga from the hazard curve
for i = 2:length(PGA)
    PPGA(i) = (1-exp(-ARE(i-1)))-(1-exp(-ARE(i)));
end
PPGA = PPGA';

r = 0.05; % discount rate
N = 50; % control period in years

%%Hospital Building
[HB_Designs,~] = xlsread('HB_Designs_sorted_short.xlsx');
HB_IC = xlsread('HB_Initial_Cost_short.xlsx');
HB_Initial_Cost = HB_IC(x(1)); % 2016 Cost Estimate

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Expected Failure Cost
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Upload1 = importdata('1Drift.txt');
Upload2 = importdata('1Damage_ratio.txt');
Drift_vs_Damage =  [Upload1.data Upload2.data];

Upload1 = importdata('HB_2Damage_ratio.txt');
Upload2 = importdata('HB_2Recovery_time.txt');
Damage_vs_Recovery =  [Upload1.data Upload2.data];

Upload1 = importdata('3Damage_ratio.txt');
Upload2 = importdata('3Minor_injury_rate.txt');
Damage_vs_Minor_injury =  [Upload1.data Upload2.data];

Upload1 = importdata('4Damage_ratio.txt');
Upload2 = importdata('4Major_injury_rate.txt');
Damage_vs_Major_injury =  [Upload1.data Upload2.data];

Upload1 = importdata('5Damage_ratio.txt');
Upload2 = importdata('5Fatality_rate.txt');
Damage_vs_Fatality =  [Upload1.data Upload2.data];

Upload1 = importdata('HB_PWMinor_drift.txt');
Upload2 = importdata('HB_PWMinor_prob.txt');
Drift_vs_PWMinor_Damage =  [Upload1.data Upload2.data];

Upload1 = importdata('HB_PWModerate_drift.txt');
Upload2 = importdata('HB_PWModerate_prob.txt');
Drift_vs_PWModerate_Damage =  [Upload1.data Upload2.data];

Upload1 = importdata('HB_PWExtensive_drift.txt');
Upload2 = importdata('HB_PWExtensive_prob.txt');
Drift_vs_PWExtensive_Damage =  [Upload1.data Upload2.data];

Upload1 = importdata('HB_PWComplete_drift.txt');
Upload2 = importdata('HB_PWComplete_prob.txt');
Drift_vs_PWComplete_Damage =  [Upload1.data Upload2.data];

% fraction converted from one level of injury to other
Upload1 = importdata('11Damage_ratio.txt');
Upload2 = importdata('11Min_Maj_transfer.txt');
Damage_vs_Min_to_Maj_injury =  [Upload1.data Upload2.data];

Upload1 = importdata('12Damage_ratio.txt');
Upload2 = importdata('12Maj_Fat_transfer.txt');
Damage_vs_Maj_to_Fatality =  [Upload1.data Upload2.data];

% Ref. MCEER report (System performance under multi-hazard environments)
DrR = [0.7  2.5  5]; % Drift Ratio limits based on FEMA 356-2000

% Amplified factor=ratio of area of structure/MCEER demonstration hospital
% structural costs
Str_RC = (138^2)/(275*56.5)*[280000 1512000 67500000  67500000]; % repair cost
Str_IC = (138^2)/(275*56.5)*[0      4380000 203670000 203670000]; % Indirect cost

% HVAC 
HVAC_Chillers = ceil((138^2)/(275*56.5)*2); % rounded up to next integer
AL_HVAC = [2 4]; % acceleration limits
HVAC_RC = [90000    500000]; % HVAC repair cost per chiller
HVAC_IC = [139500   1395000]; % HVAC indirect cost per chiller

% Partition walls
PW_RC = [230  460  690  920]; % partition wall repair cost
PW_IC = [1500 3000 4500 4500]; % partition wall indirect cost 

% Piping
DL_P =  [1.1 2.2]; % drift ratio limits for existing system
P_RC =  (138^2)/(275*56.5)*[5380 12770 14050];     % piping repair cost for all the floors combinedly
P_IC =  (138^2)/(275*56.5)*[0 97650 1046250];      % piping indirect cost all floors

%
Max_CL = 1.85*HB_Initial_Cost;   % Loss of Contents for 100% damage (1.85 factor from ATC-13, Table 4.11)
Max_RC = 276544.5573;   % Relocation Cost per month for 100% damage

Max_MinIL = 854599.5;     % Minor Injury Loss for 100% damage
Max_MajIL = 8545995;    % Major Injury Loss for 100% damage
Max_DL = 4021781613;   % Death Loss for 100% damage

% Drift Prediction
for j = 1:length(PGA)
    NN_Input(j,:) =  [sec_prop(HB_Designs(x(1),1),10) sec_prop(HB_Designs(x(1),2),10) sec_prop(HB_Designs(x(1),3),10) sec_prop(HB_Designs(x(1),4),10) ...
        sec_prop(HB_Designs(x(1),5),10) sec_prop(HB_Designs(x(1),6),10) sec_prop(HB_Designs(x(1),7),10) sec_prop(HB_Designs(x(1),8),10) 50 PGA(j)];
end 
Drift_Output = HB_NNFcn_BR_Drift_Ratio(NN_Input');
Drift_ratio_HB = round(Drift_Output',2);   

Acc_Output = HB_NNFcn_BR_Abs_Acc(NN_Input');
Abs_Acc = round(Acc_Output',2);

for j = 1:length(Drift_ratio_HB)
    if Drift_ratio_HB(j)<0.2
        DmR = 0;
        HB_RT(j) = 0;
        MinIF = 0;
        MajIF = 0;
        DF = 0;
    elseif (Drift_ratio_HB(j)>=0.2) && (Drift_ratio_HB(j)<=5)
        DmR = Drift_vs_Damage(find(Drift_ratio_HB(j)==(Drift_vs_Damage(:,1)),1),2);
        HB_RT(j) = Damage_vs_Recovery(find(DmR==(Damage_vs_Recovery(:,1)),1),2);

          initial_MinIF = Damage_vs_Minor_injury(find(DmR==(Damage_vs_Minor_injury(:,1)),1),2);
          if DmR<0.006
              Min_to_MajIF = 0;
          else
              Min_to_MajIF = Damage_vs_Min_to_Maj_injury(find(DmR==(Damage_vs_Min_to_Maj_injury(:,1)),1),2);
          end
          MinIF = initial_MinIF*(1-Min_to_MajIF/100);

          if DmR<=0.8

          initial_MajIF = Damage_vs_Major_injury(find(DmR==(Damage_vs_Major_injury(:,1)),1),2);
          if DmR<0.005
              Maj_to_FatF = 0;
          else
              Maj_to_FatF = Damage_vs_Maj_to_Fatality(find(DmR==(Damage_vs_Maj_to_Fatality(:,1)),1),2);
          end
          MajIF = initial_MajIF*(1-Maj_to_FatF/100) + initial_MinIF*Min_to_MajIF/100; 

          DF = Damage_vs_Fatality(find(DmR==(Damage_vs_Fatality(:,1)),1),2);
          DF = DF + initial_MajIF*Maj_to_FatF/100;

          else
              % linear interpolation
              initial_MajIF = 0.04+((DmR-0.8)/0.2)*(0.4-0.04);
              Maj_to_FatF = 100;
              MajIF = initial_MajIF*(1-Maj_to_FatF/100) + initial_MinIF*Min_to_MajIF/100; 

              DF = 0.01+((DmR-0.8)/0.2)*(0.2-0.01);
              DF = DF + initial_MajIF*Maj_to_FatF/100;
          end
      else
          DmR = 1;
          HB_RT(j) = 723.4;
  %       Minor injury completely transferred into major injury and major injury into death        
  %       MinIF = 0.4; 
          MinIF = 0;
          MajIF = 0.4;
          DF = 0.2+0.4;           
      end

          CL(j) = Max_CL*DmR;
          RC(j) = Max_RC*HB_RT(j)/30; 

          MinIL(j) = Max_MinIL*MinIF;
          MajIL(j) = Max_MajIL*MajIF;
          DL(j) = Max_DL*DF;

      % Structural cost evaluation
      % Performance level identification
      c = find(Drift_ratio_HB(j)<DrR,1); % to find the index of the matrix DrR
      if isempty(c)==0          
          PL = c;   % Performance Level
      else PL = 4;  % Drift > 5
      end
      Src(j) = Str_RC(PL); % repair cost
      Sic(j) = Str_IC(PL); % indirect cost  

      % HVAC cost evaluation
      % Performance level identification
      c = find(Abs_Acc(j)<AL_HVAC,1);
      if isempty(c)==0          
          PL = c;   % Performance Level
      else PL = 3;  % Acc > 4
      end
      if PL==1
          HVACrc(j) = 0;
          HVACic(j) = 0;
      else
          HVACrc(j) = HVAC_Chillers*HVAC_RC(PL-1);
          HVACic(j) = HVAC_Chillers*HVAC_IC(PL-1);
      end

      % Partition wall cost evaluation
      PWrc(j) = 0;
      PWic(j) = 0;
      beds = ceil(138^2/(275*57.5)*93);
      levels = {'Minor' 'Moderate' 'Extensive' 'Complete'};
      for l=1:length(levels)
          upload1 = importdata(char(strcat('HB_PW',levels(l),'_drift.txt')));
          upload2 = importdata(char(strcat('HB_PW',levels(l),'_prob.txt'))); 
          e = upload1.data;
          f = upload2.data; 
          g = Drift_ratio_HB(j);
          % probability of exceedance damage state
          if g<=min(e)
              pds(l) = 0;
          elseif g>max(e)
              pds(l) = 1;
          else pds(l) = f(find(g==e,1));
          end   
      end
      % probability of each damage state
      pd = pds;
      for m=1:(length(levels)-1)
          pd(m) = pds(m) - pds(m+1);
      end
      PWrc(j) = pd*PW_RC';            
      PWic(j) = beds*pd*PW_IC';

      % Piping cost evaluation
      % Performance level identification
      Prc(j)=0;
      Pic(j)=0;
      c=find(Drift_ratio_HB(j)<DL_P,1);
      if isempty(c)==0          
          PL = c;   % Performance Level
      else PL = 3;  % Drift ratio > 2.2
      end
      Prc(j) = P_RC(PL);
      Pic(j) = P_IC(PL);                   
  end

HB_Annual_DD_Cost =  (Src+HVACrc+PWrc+Prc+CL+Sic+HVACic+PWic+Pic+RC+MinIL+MajIL+DL)*PPGA; % Direct+Indirect+Social including death loss damage cost
HB_Exp_DD_Cost = HB_Annual_DD_Cost/r*(1-(1/(1+r)^N));
% HB_Rest_Time = RT*PPGA;

y(1) = HB_Initial_Cost;
y(2) = HB_Exp_DD_Cost;

[Cesar Antonio Lopez Segura]

Hi Vamshi,

To improve the efficiency of your code the first step could be separate the data processing (importdata, variable creation,... ) outside of function. And create a structure with all your data ( see the example below ):

Mydata.PGA = PGA';
Mydata. Are = Are;
Mydata.r = 0.05;
Mydata.N = 50;

Then, you can attach the structure variable "Mydata" and create a simplify function.

To use Mydata structure in your function it is possible to create a new input or declare Mydata as global.

1.-

    function y = HB_RegionalLoss_InitialCost(x, Mydata)

2.-

function y = HB_RegionalLoss_InitialCost(x)
global Mydata

Finally, you can share (attach) your function and Mydata structure with MATLAB community.

Note: it is necessary that objective function has a single output. The code example has two outputs:

% HB_Rest_Time = RT*PPGA;
y(1) = HB_Initial_Cost;
y(2) = HB_Exp_DD_Cost;

If you need to minimize two objectives the function output could be:

y = WeightOne* HB_Initial_Cost +  WeightTwo*HB_Exp_DD_Cost

The value of WeightOne and WeightTwo could depend of your priorities.

[Vamshi Krishna Gudipati]

Hello Cesar,

Thank you for your suggestions. It helped a lot to improve the run time.

[Cesar Antonio Lopez Segura]

Hi Vamshi,

Do you have the global minimum of your function ?

[Vamshi Krishna Gudipati]

Yes