how to incorporate boundary conditions as constraints in fmincon optimization?

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)

[Torsten]

I don't see where you try to define constraints.

These two lines

norm_cons = sqrt(probability2);
variational_pi1DB_fun = t_pi1DB/norm_cons;

are nowhere used in the problem definition, as far as I can see.

[Saurav Parmar]

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]

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.

Answer 1

[Torsten]

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]

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]

Maybe the documents and examples available here:

https://uk.mathworks.com/products/optimization.html

https://uk.mathworks.com/help/optim/solver-based-approach.html

https://uk.mathworks.com/help/optim/problem-based-approach.html

https://uk.mathworks.com/help/optim/optimization-problem-setup-solver-based.html

https://uk.mathworks.com/help/optim/optimization-problem-setup-problem-based.html

?

Answer 2

[Matt J]

Use fcn2optimexpr to implement nonlinear constraints.

[Saurav Parmar]

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.