how to find the values that minimize the objective function using the generic algorithm?
이전 댓글 표시
i have written the following script that finds the value of the objective function between experimental and modelling data. I want to minimize the objective function by using the generic algorithm. but it always stucks on the bounds that i set. Any suggestions on how to fix this?
clear all
close all
clc
% Define the experimental data
exp_t = [0, 5, 15, 30, 45, 60]; % time in mins
exp_y = [100, 0, 0, 0, 0, 0; % Cop yield
0, 32.3, 31.9, 31, 29.7, 28; % Bio-oil yield
0, 33.3, 33, 29.9, 23.5, 13.8; % Biochar yield
0, 24.8, 25.8, 26.8, 27.1, 26.7; % AP yield
0, 9.7, 9.2, 12.2, 19.6, 31.5]; % Gas yield
% Define the constant temperature
T = 300 + 273.15;
R = 8.314;
lb = [10 2];
ub = [100 25000];
guess = [8 10000 8 10000 8 10000 8 10000 8 10000 8 10000 8 10000];
% Preallocate k
k = zeros(1, numel(guess)/2);
for i = 1:2:numel(guess)
k((i+1)/2) = guess(i) * exp(-(guess(i+1)/(R*T)));
end
% Initialize concentrations
x0= [100, 0, 0, 0, 0]; % Initial concentrations
ode = @(t, x) [-(k(1) + k(2) + k(3) + k(4)) * x(1); % dMcop/dt
k(1) * x(1) - (k(5) + k(6)) * x(2); % dMbiooil/dt
k(2) * x(1); % dMBiochar/dt
k(3) * x(1) + k(5) * x(2) - k(7) * x(4); % dMaqueous/dt
k(4) * x(1) + k(6) * x(2) + k(7) * x(4)]; % dMgas/dt
[t, x] = ode45(ode, exp_t, x0); % solve numerically
% Calculate the squared differences
sq_diff = (exp_y - x').^2;
% Sum of squared differences
objective_function = sum(sq_diff(:));
disp(['Objective Function Value: ' num2str(objective_function)]);
% Define the objective function as a function of A and E
objfun = @(AE) objfun_internal(AE, exp_t, exp_y, T, R);
% Define the anonymous function that calls objfun
objfun2 = @(AE) objfun_internal([AE(1) AE(2) 8 10000 8 10000 8 10000 8 10000 8 10000 8 10000 8 10000], exp_t, exp_y, T, R);
% Define the internal objfun function
function result = objfun_internal(AE, exp_t, exp_y, T, R)
% Place your existing objfun logic here
k = zeros(1, numel(AE)/2);
for i = 1:2:numel(AE)
k((i+1)/2) = AE(i) * exp(-(AE(i+1)/(R*T)));
end
% Initialize concentrations
x0 = [100, 0, 0, 0, 0]; % Initial concentrations
ode = @(t, x) [-(k(1) + k(2) + k(3) + k(4)) * x(1); % dMcop/dt
k(1) * x(1) - (k(5) + k(6)) * x(2); % dMbiooil/dt
k(2) * x(1); % dMBiochar/dt
k(3) * x(1) + k(5) * x(2) - k(7) * x(4); % dMaqueous/dt
k(4) * x(1) + k(6) * x(2) + k(7) * x(4)]; % dMgas/dt
[t, x] = ode45(ode, exp_t, x0); % solve numerically
% Calculate the squared differences
sq_diff = (exp_y - x').^2;
% Sum of squared differences
result = sum(sq_diff(:));
end
댓글 수: 6
Sam Chak
2024년 1월 12일
Ciao Georgios, I don't find the genetic algorithm solver ga() in your code. Can you include it?
Hi @GEORGIOS
I see that you're attempting to determine the parameters for the biochemical system of differential equations. I've conducted my run, and the objective value remains relatively high at 13K. Is it possible that there are dynamics not accounted for in your suggested ordinary differential equations (ODEs)?
% Define the experimental data
exp_t = [0, 5, 15, 30, 45, 60]; % time in mins
exp_y = [1, 0, 0, 0, 0, 0; % Cop yield
0, 32.3, 31.9, 31, 29.7, 28; % Bio-oil yield
0, 33.3, 33, 29.9, 23.5, 13.8; % Biochar yield
0, 24.8, 25.8, 26.8, 27.1, 26.7; % AP yield
0, 9.7, 9.2, 12.2, 19.6, 31.5]; % Gas yield
% Define the constant temperature
T = 300 + 273.15;
R = 8.314;
guess = [50 100 50 100 50 100 50 100 50 100 50 100 50 100];
% Define the optimization options
options = optimoptions('ga');
% options.Display = 'iter'; % Display optimization information during iterations
options.MaxGenerations = 500; % Maximum number of generations
options.PopulationSize = 50; % Population size
options.ConstraintTolerance = 0; % Adjust constraint tolerance
% Lower bounds for A and E
lb = [1 5 1 5 1 5 1 5 1 5 1 5 1 5];
% Perform optimization using genetic algorithm
[optimal_params, min_objective] = ga(@(params) objective_function_wrapper(params, exp_t, exp_y, R, T), numel(guess), [], [], [], [], lb, [], [], options)
% Extract optimized A and E values
optimized_A = optimal_params(1:2:end);
optimized_E = optimal_params(2:2:end);
% Display results
disp(['Optimal A values: ' num2str(optimized_A)]);
disp(['Optimal E values: ' num2str(optimized_E)]);
disp(['Minimum Objective Function Value: ' num2str(min_objective)]);
% Function wrapper for the objective function
function obj = objective_function_wrapper(params, exp_t, exp_y, R, T)
k = zeros(1, numel(params)/2);
for i = 1:2:numel(params)
k((i+1)/2) = params(i) * exp(-(params(i+1)/(R*T)));
end
x0 = [1, 0, 0, 0, 0]; % Initial concentrations
ode = @(t, x) [-(k(1) + k(2) + k(3) + k(4)) * x(1);
k(1) * x(1) - (k(5) + k(6)) * x(2);
k(2) * x(1);
k(3) * x(1) + k(5) * x(2) - k(7) * x(4);
k(4) * x(1) + k(6) * x(2) + k(7) * x(4)];
[t, x] = ode45(ode, exp_t, x0);
sq_diff = (exp_y - x').^2;
obj = sum(sq_diff(:));
end
GEORGIOS
2024년 1월 12일
this schematic represents the process that happens, assuming all reactions are first order:
the data about the yields that i have are these:
I have adjusted them and added an extra line for t=0 as initial conditions. I am not exactly sure if this reasoning was correct though.
Your measurement data are in wt %. So there are two things you should ask yourself:
- Usually, kinetic models work in mol %, not in wt %. Is this different in your case ?
- Why the sums of each row do not add to 100 ?
And a final remark:
My guess is that you will have to divide your experimental data by a factor of 100 (except for the first row). Your model works in fractions, not in percent.
GEORGIOS
2024년 1월 13일
답변 (2개)
Hassaan
2024년 1월 12일
0 개 추천
- Adjust Bounds: Review the lower (lb) and upper (ub) bounds that you've set for the optimization variables (A and E). Make sure that the bounds are reasonable and wide enough to allow the optimization algorithm to explore the search space effectively. Tight bounds can cause the optimization to get stuck.
- Population Size and Generations: Experiment with the population size and the number of generations in the genetic algorithm. Increasing the population size and/or the number of generations can help the algorithm explore the search space more thoroughly.
- Mutation and Crossover Rates: Adjust the mutation rate and crossover rate parameters of the genetic algorithm. These parameters control the exploration-exploitation trade-off. Higher mutation rates can help escape local minima, while lower mutation rates can help convergence. Experiment with different values to find a balance.
- Objective Function Scaling: Consider scaling your objective function. Sometimes, optimization algorithms work better when the objective function values are within a certain range. Scaling can help in achieving this.
- Initial Guess: Experiment with different initial guesses for the optimization variables. The algorithm's convergence can depend on the starting point, and different initial points can lead to different local minima.
- Convergence Criteria: Check the convergence criteria of the genetic algorithm. Ensure that the algorithm is set to terminate properly when it converges or when a maximum number of function evaluations is reached.
- Parallel Computing: If you have access to multiple CPU cores, consider enabling parallel computing for the genetic algorithm. This can speed up the optimization process.
- Optimization Toolbox Functions: MATLAB's Optimization Toolbox provides various optimization functions that are tailored for different optimization problems. You might explore other optimization techniques such as fmincon, fminunc, or patternsearch in addition to the genetic algorithm.
- Output Verbosity: Increase the verbosity of the optimization output. This can provide insights into what's happening during optimization and help you identify potential issues.
- Plot the Optimization Progress: Visualize the optimization progress by plotting the objective function values over iterations or generations. This can help you understand the optimization dynamics.
- Advanced Techniques: If the optimization problem is complex and standard optimization methods do not work well, you might explore more advanced techniques like surrogate modeling or global optimization algorithms.
---------------------------------------------------------------------------------------------------------------------------------------------------------
If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.
Professional Interests
- Technical Services and Consulting
- Embedded Systems | Firmware Developement | Simulations
- Electrical and Electronics Engineering
Feel free to contact me.
Star Strider
2024년 1월 12일
I used some standardised code that I wrote for these sorts of problems that I adapted to your data —
clear all
close all
clc
% Define the experimental data
exp_t = [0, 5, 15, 30, 45, 60]; % time in mins
exp_y = [1, 0, 0, 0, 0, 0; % Cop yield
0, 32.3, 31.9, 31, 29.7, 28; % Bio-oil yield
0, 33.3, 33, 29.9, 23.5, 13.8; % Biochar yield
0, 24.8, 25.8, 26.8, 27.1, 26.7; % AP yield
0, 9.7, 9.2, 12.2, 19.6, 31.5]; % Gas yield
t = exp_t.';
c = exp_y.';
% % Define the constant temperature
% T = 300 + 273.15;
% R = 8.314;
% lb = [10 2];
% ub = [100 25000];
% guess = [8 10000 8 10000 8 10000 8 10000 8 10000 8 10000 8 10000];
% % Preallocate k
% k = zeros(1, numel(guess)/2);
% for i = 1:2:numel(guess)
% k((i+1)/2) = guess(i) * exp(-(guess(i+1)/(R*T)));
% end
% % Initialize concentrations
% x0= [1, 0, 0, 0, 0]; % Initial concentrations
% ode = @(t, x) [-(k(1) + k(2) + k(3) + k(4)) * x(1); % dMcop/dt
% k(1) * x(1) - (k(5) + k(6)) * x(2); % dMbiooil/dt
% k(2) * x(1); % dMBiochar/dt
% k(3) * x(1) + k(5) * x(2) - k(7) * x(4); % dMaqueous/dt
% k(4) * x(1) + k(6) * x(2) + k(7) * x(4)]; % dMgas/dt
% [t, x] = ode45(ode, exp_t, x0); % solve numerically
% % Calculate the squared differences
% sq_diff = (exp_y - x').^2;
% % Sum of squared differences
% objective_function = sum(sq_diff(:));
% disp(['Objective Function Value: ' num2str(objective_function)]);
% % Define the objective function as a function of A and E
% objfun = @(AE) objfun_internal(AE, exp_t, exp_y, T, R);
% % Define the anonymous function that calls objfun
% objfun2 = @(AE) objfun_internal([AE(1) AE(2) 8 10000 8 10000 8 10000 8 10000 8 10000 8 10000 8 10000], exp_t, exp_y, T, R);
% % Define the internal objfun function
% function result = objfun_internal(AE, exp_t, exp_y, T, R)
% % Place your existing objfun logic here
% k = zeros(1, numel(AE)/2);
% for i = 1:2:numel(AE)
% k((i+1)/2) = AE(i) * exp(-(AE(i+1)/(R*T)));
% end
% % Initialize concentrations
% x0 = [1, 0, 0, 0, 0]; % Initial concentrations
% ode = @(t, x) [-(k(1) + k(2) + k(3) + k(4)) * x(1); % dMcop/dt
% k(1) * x(1) - (k(5) + k(6)) * x(2); % dMbiooil/dt
% k(2) * x(1); % dMBiochar/dt
% k(3) * x(1) + k(5) * x(2) - k(7) * x(4); % dMaqueous/dt
% k(4) * x(1) + k(6) * x(2) + k(7) * x(4)]; % dMgas/dt
% [t, x] = ode45(ode, exp_t, x0); % solve numerically
% % Calculate the squared differences
% sq_diff = (exp_y - x').^2;
% % Sum of squared differences
% result = sum(sq_diff(:));
% end
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 500;
Parms = 12;
optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10); % Options Structure For 'Answers' Problems
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
optshf = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'HybridFcn',@fmincon, 'PlotFcn',@gaplotbestf, 'PlotInterval',1); % With 'HybridFcn'
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts);
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)
fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','N')
function C=kinetics(k,t)
c0 = k(8:12);
% c0=[1;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,x)
dC = [-(k(1) + k(2) + k(3) + k(4)) * x(1); % dMcop/dt
k(1) * x(1) - (k(5) + k(6)) * x(2); % dMbiooil/dt
k(2) * x(1); % dMBiochar/dt
k(3) * x(1) + k(5) * x(2) - k(7) * x(4); % dMaqueous/dt
k(4) * x(1) + k(6) * x(2) + k(7) * x(4)]; % dMgas/dt
end
C=Cv;
end
This requires more than the 55 seconds allowed here, so you will need to run it offline.
This code runs and gives reasonable result, however it does not give a good fit to the data, even on several runs. It may be necessary to review your system of differential equations to be sure that they correctly describe the problem you want to solve. (With other data and other problems, my code here generally gives quite good results.) Note that it also estimates the initial conditions of the differential equations, and reports those as the last 5 elements of the 12-elrment ‘theta’ parameter vector.
I do not completely udnerstand what you are doing, so it will be necessary for you to add the calculations and such involving temperature to get the result you want.
I have this running in the background (so far for 20 minutes). When it converges, I will post its results.
.
댓글 수: 7
GEORGIOS
2024년 1월 12일
thanks for your answer. I cannot exactly understand what your code is doing. I want to calculate the activation energies and the pre exponential factor of the Arhennius equation, which you completely commented out, and that then will be used to solve the equations. Furthermore, do you use the temperature value somewhere in this code? Last, i do not believe that this problem is so complex to run for hours since its for a masters and the proffessor hasnt mentioned something about this being so complex and that would require so much time to run.
The following picture shows a good representation of what i want my code to do:
Star Strider
2024년 1월 12일
My pleasure.
My objective was to get the ga part of the code working. My code uses ga rather than any genetic algorithm code I wrote.
The essence of it is to initialise the fitness function and the initial population matrix:
t = exp_t;
c = exp_y;
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 500;
Parms = 12;
then call ga:
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts);
however first set up the objective function ‘kinetics’ using the integrated differential equations:
function C=kinetics(k,t)
c0 = k(8:12);
% c0=[1;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,x)
dC = [-(k(1) + k(2) + k(3) + k(4)) * x(1); % dMcop/dt
k(1) * x(1) - (k(5) + k(6)) * x(2); % dMbiooil/dt
k(2) * x(1); % dMBiochar/dt
k(3) * x(1) + k(5) * x(2) - k(7) * x(4); % dMaqueous/dt
k(4) * x(1) + k(6) * x(2) + k(7) * x(4)]; % dMgas/dt
end
C=Cv;
end
The rest of it just reports and plots the data.
I did not use any of the conversions because I do not understand what you are doing. If those need to be part of the objective function (‘kinetics’), then add them as appropriate to it so that it reports out the correct matrix as ‘C’. If you want to add temperature and anything else to ‘kinetics’ declare them as additional arguments:
function C=kinetics(theta,t,T)
...
end
Since ‘C’ (the output of ‘kinetics’) is what gets compared to your data, it needs to match as closely as possible to them. If that requires additional calculations, do them in ‘kinetics’.
The rest of the code does not need to be changed even if ‘kinetics’ is changed internally to add more calculations to the result.
After letting it run for nearly 40 minutes, it is still not converging on the correct solution. It stopped when it thought it had converged.
GEORGIOS
2024년 1월 21일
i also wrote this simpler code but the fit isn't good:
clear all
close all
clc
% Define the experimental data for different temperatures
t = [0, 5, 15, 30, 45, 60]; % time in mins
% Example: Create cell arrays for exp_y and T for three different temperatures
exp_y= [1 0 0 0 0;
0 0.299 0.361 0.21 0.13;
0 0.299 0.344 0.224 0.133;
0 0.295 0.289 0.239 0.177;
0 0.288 0.202 0.247 0.263;
0 0.276 0.082 0.248 0.394];
c1 = exp_y(:,1);
c1_0 = 1;
c2 = exp_y(:,2);
c2_0 = 0;
c3 = exp_y(:,3);
c3_0 = 0;
c4 = exp_y(:,4);
c4_0 = 0;
c5 = exp_y(:,5);
c5_0 = 0;
tin = t';
yin = [c1,c2,c3,c4,c5];
p0 = [0.027 0.0016 0.00023 0.0013 0.004 0.013 0.061 0.080 c1(1) c2(1) c3(1) c4(1) c5(1)]; % initial guess
lb = [0 0 0 0 0 0 0 0 0 0 0 0 0];
ub = [0.48 0.48 0.48 0.48 0.48 0.48 0.48 0.48 1 1 1 1 1];
%Lsqcurvefit
options = optimoptions('lsqcurvefit','StepTolerance',1e-100,'MaxFunctionEvaluations',10000,'MaxIterations',10000);
[p,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@(p,time) diff1(p,time,yin),p0,tin,yin,lb,ub,options);
FitData = diff1(p,t);
c1_out = FitData(:,1);
c2_out = FitData(:,2);
c3_out = FitData(:,3);
c4_out = FitData(:,4);
c5_out = FitData(:,5);
time_plot = [0,t];
c1_plot = [c1_0,c1'];
c2_plot = [c2_0,c2'];
c3_plot = [c3_0,c3'];
c4_plot = [c4_0,c4'];
c5_plot = [c5_0,c5'];
%Plot
c1_out_plot = [c1_0,c1_out'];
c2_out_plot = [c2_0,c2_out'];
c3_out_plot = [c3_0,c3_out'];
figure(2)
plot(time_plot,c1_plot,'bo')
grid on
hold on
plot(time_plot,c2_plot,'ro')
plot(time_plot,c3_plot, 'go')
plot(time_plot,c1_out_plot,'b')
plot(time_plot,c2_out_plot,'r')
plot(time_plot,c3_out_plot,'g')
hold off
xlabel('Time [min]')
ylabel('Concentration (wt.%)')
legend('c1','c2', 'c3')
SStotC1 = sum((c1 - mean(c1)).^2);
SSresC1 = sum((c1(:) - c1_out).^2);
R2C1 = 1 - SSresC1/SStotC1
SStotC2 = sum((c2 - mean(c2)).^2);
SSresC2 = sum((c2(:) - c2_out).^2);
R2C2 = 1 - SSresC2/SStotC2
SStotC3 = sum((c3 - mean(c3)).^2);
SSresC3 = sum((c3(:) - c3_out).^2);
R2C3 = 1 - SSresC3/SStotC3
SStotC4 = sum((c4 - mean(c4)).^2);
SSresC4 = sum((c4(:) - c4_out).^2);
R2C4 = 1 - SSresC4/SStotC4
SStotC5 = sum((c5 - mean(c5)).^2);
SSresC5 = sum((c5(:) - c5_out).^2);
R2C5 = 1 - SSresC5/SStotC5
p(1:8)
function C =diff1(p,time,~)
x0 = p(9:13);
[t,Cout] = ode45(@DifEq,time,x0);
function dC = DifEq(t,x)
xdot = zeros(5,1);
xdot(1) = -(p(1) + p(2) + p(3)) * x(1);
xdot(2) = p(2) * x(1) + p(4) * x(4) - (p(5) + p(6) + p(8)) * x(2);
xdot(3) = p(3) * x(1) + p(6) * x(2);
xdot(4) = p(1) * x(1) + p(5) * x(5) - (p(4) + p(7)) * x(4);
xdot(5) = p(8) * x(2) + p(7) * x(4);
dC = xdot;
end
C = Cout;
end
Star Strider
2024년 1월 21일
My pleasure!
I will re-run my code with your new data later.
‘... can i set different bounds for parameters in your code?’
I will look at the rest of that (especially your latest comment) later.
I do not understand what you are doing with this (and similar sets of statements):
SStotC1 = sum((c1 - mean(c1)).^2);
SSresC1 = sum((c1(:) - c1_out).^2);
R2C1 = 1 - SSresC1/SStotC1
My ‘ftns’ function can be made to calculate those if that is what you want to do.
.
As always, you have much too few data for far too many coefficients to be estimated. Mathematics cannot substitute the lack of data.
And therefore the temperature influence on your conversion is far too small to be seen in parameter trends.
If you implicitly want to implement that the rate constants increase with increasing temperature, fit all data simultaneously and use the Arrhenius law for the rate expressions which are temperature-dependent (although this again will introduce activation energies as new unknown parameters).
카테고리
도움말 센터 및 File Exchange에서 Nonlinear Optimization에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
