I have two problems in my code: 1 - Genetic algorithm gets “stuck” in the first generation. 2 – lsqcurvefit is unable to fix it. Can anyone, please, help me?

Hello, I am trying to minimize the error between 2 curves using the least squares method. To accomplish that I am using a two-step fitting procedure with the GA and the lsqcurvefit. I have some kind of issue that I cannot fix: The genetic algorithm gets stuck in the first generation, i.e. the stopping criteria is always the stall generation number, no matter the changes I make in my code. As for the lsqcurvefit, I was hoping that this procedure could go around this GA issue but I was mistaken. Probably I am badly using the gaoptimset and optimset options, but I can’t figure it out what is missing. Many thanks in advance for any help.

I have 9 variables that I need to optimize. For each variable I have gave it a lower and an upper bound.

% Lower bound

LB = horzcat(15,5,20,1E5,1,1E7,0,0,0);

% Upper bound

UB = horzcat(25,15,40,1E7,1E2,1E8,1,1,1);

I am trying in a test case and I know that the optimal values are, in the same order:

19.1, 10.0, 31.0, 2E6, 4, 6.5E7, 0, 1, 0

GA code:

delta = @FittingObj.Fitting_GA;

nvars = numel(LB);

options_GA = gaoptimset('PlotFcn',{@gaplotrange,@gaplotbestf,@gaplotdistance});
options_GA = gaoptimset(options_GA,'PopInitRange',[LB;UB]); 
options_GA = gaoptimset(options_GA,'Generations',150,'StallGenLimit',20); 
options_GA = gaoptimset(options_GA,'PopulationSize',20); 
% I tried both the 'StallGenLimit' and the 'PopulationSize' with 100 and the result is the same
options_GA = gaoptimset(options_GA,'SelectionFcn',{@selectiontournament,4});
options_GA = gaoptimset(options_GA,'CrossoverFcn',{@crossoverscattered});
options_GA = gaoptimset(options_GA,'MutationFcn',{@mutationuniform,0.10});
options_GA = gaoptimset(options_GA,'Display','iter');

% Run the GA solver.
[y,opt_delta_GA] = ga(delta,nvars,[],[],[],[],LB,UB,[],options_GA);

% And my objective function is:

function delta_GA = Fitting_GA(fit,y)            
            fE2 = y(1);
            fE3 = y(2);
            fE4 = y(3);
            fA2 = y(4);
            fA3 = y(5);
            fA4 = y(6);
            fb2 = y(7);
            fb3 = y(8);
            fb4 = y(9);

              fit.ReWriteFile(fE2,fE3,fE4,fA2,fA3,fA4,fb2,fb3,fb4);
              fit.pyr.InitPyrolysis();          
              fit.pyr.RunPyrolysis(); 
              fit.RearrangeDimensions();

              delta_2D = zeros(1,fit.pyr.NI);
              delta_2D(1:fit.pyr.NI) = (fit.wt_loss_ref(1:fit.pyr.NI)-fit.pyr.wt_loss(1:fit.pyr.NI)).^2;

              delta_1D = sum(delta_2D);
              fit.delta_wt_loss = (fit.pyr.NI^(-1)*delta_1D)^(1/2);
              delta_GA = log(1+10*fit.delta_wt_loss); 
    end

lsqcurvefit code:

% Run LSQ solver.
X0 = y(1:nvars);
options_LSQ = optimset('MaxFunEvals',1E5,'MaxIter',1E5);
[x,fit_error] = lsqcurvefit(@FittingObj.Fitting_LSQ,X0,PyrolysisObj.solver.x,FittingObj.wt_loss_ref,LB,UB,options_LSQ); 
opt_delta = sqrt(fit_error/PyrolysisObj.NI);

% And my objective function is:

          function mp = Fitting _LSQ(fit,x,~)
              fE2 = x(1);
              fE3 = x(2);
              fE4 = x(3);
              fA2 = x(4);
              fA3 = x(5);
              fA4 = x(6);
              fb2 = x(7);
              fb3 = x(8);
              fb4 = x(9);

              fit.ReWrite(fE2,fE3,fE4,fA2,fA3,fA4,fb2,fb3,fb4);
              fit.pyr.InitPyrolysis();            
              fit.pyr.RunPyrolysis();

              mp = fit.pyr.wt_loss(1:fit.pyr.NI);
          end

and I get:

Initial point is a local minimum. Optimization completed because the size of the gradient at the initial point is less than the default value of the function tolerance.

And this is both the GA and lsqcurvefit optimum result:

15, 5, 20, 1E5,1,1E7,0,0.5528,1 (the first variables set chosen by GA)

I saved all the GA data in a .txt file that is attached.

Answer 1

[Alan Weiss]

It is possible that lsqcurvefit gets stuck because you are trying to optimize a simulation: the tiny finite difference steps that lsqcurvefit makes in trying to estimate a gradient can lead to the simulation not changing its results at all. For suggestions about optimizing a simulation, see Optimizing a Simulation or ODE, especially the suggestions on finite differences.

Of course, I am just guessing, you didn't say that you were trying to optimize a simulation, it just looks that way to me.

And might I suggest that once you get lsqcurvefit working, ditch ga, and to find a global minimum, just run lsqcurvefit starting from multiple random points in the bounds, such as

 x0 = lb + rand(size(lb)).*(ub-lb);

You could even just use MultiStart to do that.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Isabel]

Thank you for your answer Alan Weiss

Yes you are right, I am trying to optimize the simulation. I should have written that in my first post. I am already using the ode15s solver inside the RunPyrolysis() function:

    function dydt = FunctionOde(pyr,t,y)
      dydt = zeros(pyr.Smax,1);
      pyr.wdot = zeros(pyr.Smax,1);
      pyr.UpdateReactor(t,y);
      pyr.wdot(1:pyr.Smax) = netProdRates(pyr.biopyrolysis);
      for k = 1:pyr.Smax
        pyr.wdot(k) = pyr.wdot(k).*pyr.Mw(k);
        dydt(k) = pyr.wdot(k);
      end
    end

    function RunPyrolysis(pyr)
      tspan = [0 pyr.cond.time_final];
      options = odeset('RelTol',1.e-5,'AbsTol',1.e-6);
      options = odeset(options,'MaxStep',pyr.cond.Max_Time_Step);
      pyr.solver = ode15s(@pyr.FunctionOde,tspan,pyr.SpecPoliMi,options);
      [~,pyr.NI] = size(pyr.solver.y);
      %...
    end

I even tried using the patternsearch procedure (with the same function that I wrote for the GA) using that same expression that you posted to determine the initial guess (I got it from one of your answers to another post actually, thank you for that as well :P), but get very tendentious results (in attachment) and a solution very far from optimum.

The optimum values of the patternsearch are: 15, 5, 20, 2.1E6, 24,25, 2.6E7, 0.0416, 0.3301, 0.0262

here is my PS code:

delta = @FittingObj.Fitting_PS;

% Lower bound
LB = horzcat(15,5,20,1E5,1,1E7,0,0,0);
% Upper bound
UB = horzcat(25,15,40,1E7,100,1E8,1,1,1);
% Number of variables
nvars = numel(LB);

% Run patternsearch solver.
opts =  psoptimset('PlotFcn',{@psplotbestf,@psplotfuncount,@psplotmeshsize});
opts = psoptimset(opts,'MaxFunEvals',1E4,'MaxIter',1E4);
% cache -> avoid polling equal set of numbers
opts = psoptimset(opts,'SearchMethod',{@positivebasisnp1},'CompletePoll','on','Cache','on');
opts = psoptimset(opts,'InitialMeshSize',5,'MaxMeshSize',10);
opts = psoptimset(opts,'TolX',1E-12,'TolBind',1E-12);

X0 = zeros(1,nvars);
for i =1:nvars
    X0(i)= LB(i) + rand(size(LB(i))).*(UB(i)-LB(i));
end

[x,delta_PA] = patternsearch(delta,X0,[],[],[],[],LB,UB,opts);

It just doesn't get there... I don't understand why.

[Alan Weiss]

I can see several things in your code that are suboptimal, but do not know why your problem causes convergence to such poor values.

I suggest that you scale your problem better. x(4) should have bounds from 1 to 100, and internally you should take 1e5*x(4). Similarly, x(6) should go from 1 to 10, and internally take 1e7*x(6). In this way, all of your problem parameters are from 0 through 100, and all ranges of values are from 1 to 100.

If you are using patternsearch, I suggest that you not use the cache option. And you are not taking the start points correctly in patternsearch. The call should be something like this:

 % you will have to input an appropriate n
 x = zeros(n,nvars);delta_PA = zeros(n,1);
 for i = 1:n
    X0 = LB + rand(size(LB)).*(UB-LB);
    [x(i,:),delta_PA(i)] = patternsearch(delta,X0,[],[],[],[],LB,UB,[],opts);
 end

You should not set the ode15s options inside the function call. Set the options once in a variable, then pass the variable in, something like

options = odeset('RelTol',1.e-5,'AbsTol',1.e-6);
options = odeset(options,'MaxStep',pyr.cond.Max_Time_Step);

function RunPyrolysis(pyr,options)
...
end

Call RunPyrolysis as @(pyr)RunPyrolysis(pyr,options). See Passing Extra Parameters.

I think that you have set your MaxFunEvals and MaxIter options too low for a 9-variable problem.

But mainly I think that you should use lsqcurvefit with appropriate finite differences. That would almost certainly be the most efficient way. You can start lsqcurvefit from a variety of initial points, as I just showed you to do with patternsearch.

By the way, did you try starting patternsearch from the point that you claim is best? If so, what happens?

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Isabel]

Thank you again for your time and answer Alen Weiss. I've located the bug but I still can't fix it.

I simplified the problem to just 1 parameter (fE2). This parameter function is to rewrite a .xml file with a mechanism that needs to use the cantera's library.

I know that the .xml file is updating, but the the ODE solver gives the exact same solution every time, ignoring the new information.

Here is a code update:

function dydt = FunctionOde(pyr,t,y)

              pyr.biopyrolysis  = IdealGasMix('biomass_K.xml');
              pyr.p = pressure(pyr.biopyrolysis);
              pyr.Mw(1:pyr.Smax) = molarMasses(pyr.biopyrolysis);
              T = interp1(pyr.cond.time_exp,pyr.cond.T_exp,t);
              pyr.composition = sprintf(pyr.formatSpec,y);
              set(pyr.biopyrolysis,'T',T,'P',pyr.p,'Y',pyr.composition);

    % here is where I print the y values to check if they are updating

              dydt = zeros(pyr.Smax,1);
              pyr.wdot = zeros(pyr.Smax,1);

              pyr.wdot(1:pyr.Smax) = netProdRates(pyr.biopyrolysis);
              for k = 1:pyr.Smax
                  pyr.wdot(k) = pyr.wdot(k).*pyr.Mw(k);
                  dydt(k) = pyr.wdot(k);
              end
          end

          function RunPyrolysis(pyr,fE2)

              pyr.Ea2 = num2str(fE2);  

              delete('Biomass_K.xml');
              xDoc = xmlread(fullfile(('Biomass_K_unchanged.xml')));

              % items counting starts on 0, always change n-1 item
              pyr.all_E_Items = xDoc.getElementsByTagName('E');
              pyr.E2_Item=pyr.all_E_Items.item(1);
              pyr.E2_Item.getFirstChild.setData(pyr.Ea2);            
              xmlwrite('Biomass_K.xml',xDoc);   

              tspan = [0 pyr.cond.time_final];
              options = odeset('RelTol',1.e-5,'AbsTol',1.e-6);
              options = odeset(options,'MaxStep',pyr.cond.Max_Time_Step);
              pyr.solver = ode15s(@pyr.FunctionOde,tspan,pyr.SpecPoliMi,options);
              [~,pyr.NI] = size(pyr.solver.y);
      end

Thank you for all the help!

[Alan Weiss]

Sorry, I don't know how to help you any further. It sounds as if you have diagnosed the problem, but don't know how to fix it. I suggest that you create another question for the forum on just this single issue.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Isabel]

Thank you for everything :)

[Isabel]

To whom it may interest, I managed to solve the problem with the cleanup() function.

It is my first line of code in the ODE function:

function dydt = FunctionOde(pyr,t,y)
  cleanup();
  pyr.biopyrolysis  = IdealGasMix('biomass_K.xml');
  %...
end