# how to incorporate boundary conditions as constraints in fmincon optimization?

7 views (last 30 days)
Saurav Parmar on 7 Dec 2023
Edited: Torsten on 7 Dec 2023
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 CommentsShow 1 older commentHide 1 older comment
Saurav Parmar on 7 Dec 2023
I tried defining constraints in % Constraints block ...it is just below the % opt variables block(top). Everything after %New varables is for the definition of Optim.problem..
Torsten on 7 Dec 2023
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 on 7 Dec 2023
Edited: Torsten on 7 Dec 2023
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))
Local minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the value of the optimality tolerance, and constraints are satisfied to within the value of the constraint tolerance.
sol = 6×1
-0.0000 -24.1985 -0.6429 2.1927 -0.2191 -0.0000
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
Saurav Parmar on 7 Dec 2023
Thanks for your insightful answer...Can I ask you for some solid resources for MATLAB optimizaton routine learning/training. It looks like I am somewhat good at basic MATLAB but woefully unprepared for Optimization problems. Thanks again!!
Torsten on 7 Dec 2023
Edited: Torsten on 7 Dec 2023

Matt J on 7 Dec 2023
Use fcn2optimexpr to implement nonlinear constraints.
Saurav Parmar on 7 Dec 2023
Thanks Matt for this option. It led me know that some variables had to be defined as symbolic in my problem. Ans I was completely unaware of this function. Thanks.

### Categories

Find more on Systems of Nonlinear Equations in Help Center and File Exchange

R2021b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!