how to incorporate boundary conditions as constraints in fmincon optimization?
7 次查看(过去 30 天)
显示 更早的评论
Trying to optimize the six parameters in a trial non-linear wavefunction for particle in a box to get the minimum energy. I trying to use fmincon. but i am getting an error that says - "Unrecognized function or variable 'variational_pi1DB_fun'." Somehow I am not able to incorporate non-linear constarints in the optimization problem. Can someone help me please! Pleas let me know if other methods of optimization is more helpful!
clear
clc
%defining optimization variables and an optimization problem object.
a = optimvar('a','LowerBound',30,"UpperBound",30);
b = optimvar('b','LowerBound',30,"UpperBound",30);
c = optimvar('c','LowerBound',30,"UpperBound",30);
d = optimvar('d','LowerBound',30,"UpperBound",30);
e = optimvar('e','LowerBound',30,"UpperBound",30);
f = optimvar('f','LowerBound',30,"UpperBound",30);
prob = optimproblem;
L = 5;
% w=1;
% v=1.5;
% constraints
% wavefunction zero @ x = 0 and x = L;
cons1 = double(subs(variational_pi1DB_fun,x,0)) == 0;
cons2 = double(subs(variational_pi1DB_fun,x,L)) == 0;
prob.Constraints.cons1 = cons1;
prob.Constraints.cons2 = cons2;
x0.a = -0.00006;
x0.b = -19;
x0.c = -9;
x0.d = 6;
x0.e = -0.90;
x0.f = 0.043;
%new variables
t_pi1DB = a + b*x + c*x^2 + d*x^3 + e*x^4 + f*x^5;
% equation 4.2 nominator
t_pi1DB_D1 = diff(t_pi1DB,x); % 1st derivative
t_pi1DB_D2 = diff(t_pi1DB_D1,x); % 2nd derivative
t_pi1DB_DD = ((-hbar^2)./(2*m)).*t_pi1DB_D2; % applying the hamiltonian to trial wfn Psi
t_pi1DB_c = (t_pi1DB)'; % conjugate of psi
t_pi1DB_prod = t_pi1DB_c.*t_pi1DB_DD; % multiplying the Psi* from left side
% equation 4.2 denominator
t_pi1DB_prod2 = t_pi1DB_c.*t_pi1DB; % hamiltonian not applied
% integration in nominator % ground state
probability1 = int(real(t_pi1DB_prod),x,0,L); % integrate from 0 to L
disp('probability = '),disp(probability1)
probability1 = subs(probability1,{hbar,m,L},{1,1,5}); % evaluate the probability at n = 1
disp('probability1 = '),disp(probability1)
% integration in denominator
probability2 = int(real(t_pi1DB_prod2),x,0,L);
probability2 = subs(probability2,L,5); % evaluate the probability at n = 1
disp('probability2 = '),disp(probability2)
norm_cons = sqrt(probability2);
variational_pi1DB_fun = t_pi1DB/norm_cons;
% application of equation 4.2
variational_Energy = probability1/probability2;
%objective function as an expression in the optimization variables.
P = variational_Energy;
%the objective function in prob.
prob.Objective = P;
show(prob)
sol = solve(prob,x0);
disp('sol = '),disp(sol)
3 个评论
Torsten
2023-12-7
But you can't use a variable before you define it.
In the constraint section
% constraints
% wavefunction zero @ x = 0 and x = L;
cons1 = double(subs(variational_pi1DB_fun,x,0)) == 0;
cons2 = double(subs(variational_pi1DB_fun,x,L)) == 0;
variational_pi1DB_fun is not yet defined.
采纳的回答
Torsten
2023-12-7
编辑:Torsten
2023-12-7
I think with the mixture of numerical, symbolic and optimization variables it is easier to set up the problem directly.
x0 = [ -0.00006; -19; -9; 6; -0.90; 0.043];
L = 5;
hbar = 1;
m = 1;
sol = fmincon(@(x)obj(x,L,hbar,m),x0,[],[],[],[],[],[],@(x)nonlcon(x,L,hbar,m))
function variational_Energy = obj(xnum,L,hbar,m)
a = xnum(1);
b = xnum(2);
c = xnum(3);
d = xnum(4);
e = xnum(5);
f = xnum(6);
syms x
%new variables
t_pi1DB = a + b*x + c*x^2 + d*x^3 + e*x^4 + f*x^5;
% equation 4.2 nominator
t_pi1DB_D1 = diff(t_pi1DB,x); % 1st derivative
t_pi1DB_D2 = diff(t_pi1DB_D1,x); % 2nd derivative
t_pi1DB_DD = ((-hbar^2)./(2*m)).*t_pi1DB_D2; % applying the hamiltonian to trial wfn Psi
t_pi1DB_c = (t_pi1DB)'; % conjugate of psi
t_pi1DB_prod = t_pi1DB_c.*t_pi1DB_DD; % multiplying the Psi* from left side
% equation 4.2 denominator
t_pi1DB_prod2 = t_pi1DB_c.*t_pi1DB; % hamiltonian not applied
% integration in nominator % ground state
probability1 = int(real(t_pi1DB_prod),x,0,L); % integrate from 0 to L
% integration in denominator
probability2 = int(real(t_pi1DB_prod2),x,0,L);
% application of equation 4.2
variational_Energy = probability1/probability2;
variational_Energy = double(variational_Energy);
end
function [C,Ceq] = nonlcon(xnum,L,hbar,m)
a = xnum(1);
b = xnum(2);
c = xnum(3);
d = xnum(4);
e = xnum(5);
f = xnum(6);
syms x
%new variables
t_pi1DB = a + b*x + c*x^2 + d*x^3 + e*x^4 + f*x^5;
t_pi1DB_c = (t_pi1DB)'; % conjugate of psi
% equation 4.2 denominator
t_pi1DB_prod2 = t_pi1DB_c.*t_pi1DB; % hamiltonian not applied
% integration in denominator
probability2 = int(real(t_pi1DB_prod2),x,0,L);
% constraints
% wavefunction zero @ x = 0 and x = L;
norm_cons = sqrt(probability2);
variational_pi1DB_fun = t_pi1DB/norm_cons;
Ceq(1) = double(subs(variational_pi1DB_fun,x,0));
Ceq(2) = double(subs(variational_pi1DB_fun,x,L));
C = [];
end
更多回答(1 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!