Alright, perfect! I took one of the examples in the Optizelle code and modified it to solve a sequence of somewhat similar problems. First, create a file called ip_resolve.m with
function ip_resolve(fname)
% Grab the Optizelle library
global Optizelle;
setupOptizelle();
% Generate an initial guess
x = [2.1;1.1];
% Allocate memory for the inequality multiplier
z = [0.;0.];
% Create an optimization state
state=Optizelle.InequalityConstrained.State.t( ...
Optizelle.Rm,Optizelle.Rm,x,z);
% Read the parameters from file
state=Optizelle.json.InequalityConstrained.read( ...
Optizelle.Rm,Optizelle.Rm,fname,state);
% Run a couple of different, similar optimization problems
for c=linspace(1,-1,10)
% Create a bundle of functions
fns=Optizelle.InequalityConstrained.Functions.t;
fns.f=MyObj();
fns.h=MyIneq(c);
% Protect against having non interior starting points
if fns.h.eval(state.x) <= 0
error('Must start optimization with a strictly interior point')
end
% Solve the optimization problem
state=Optizelle.InequalityConstrained.Algorithms.getMin( ...
Optizelle.Rm,Optizelle.Rm,Optizelle.Messaging.stdout,fns,state);
% Print the optimal solution
fprintf('The optimal point for c=%e is: (%e,%e)\n', ...
c,state.x(1),state.x(2));
% Tell the optimization that it didn't really converge
state.opt_stop = Optizelle.OptimizationStop.NotConverged;
% Reset the iteration back to 1
state.iter = 1;
% Increase the barrier parameter to something larger
state.mu = 1.;
end
end
% Squares its input
function z = sq(x)
z=x*x;
end
% Define a simple objective where
%
% f(x,y)=(x+1)^2+(y+1)^2
%
function self = MyObj()
% Evaluation
self.eval = @(x) sq(x(1)+1.)+sq(x(2)+1.);
% Gradient
self.grad = @(x) [
2.*x(1)+2.;
2.*x(2)+2.];
% Hessian-vector product
self.hessvec = @(x,dx) [
2.*dx(1);
2.*dx(2)];
end
% Define simple inequalities
%
% h(x,y)= [ x + 2y >= c ]
% [ 2x + y >= 1 ]
%
function self = MyIneq(c)
% z=h(x)
self.eval = @(x) [
x(1)+2.*x(2)-c ;
2.*x(1)+x(2)-1. ];
% z=h'(x)dx
self.p = @(x,dx) [
dx(1)+2.*dx(2) ;
2.*dx(1)+dx(2) ];
% xhat=h'(x)*dz
self.ps = @(x,dz) [
dz(1)+2.*dz(2) ;
2.*dz(1)+dz(2) ];
% xhat=(h''(x)dx)*dz
self.pps = @(x,dx,dz) zeros(2,1);
end
Then another called tr_newton.json with:
{
"Optizelle" : {
"msg_level" : 1,
"H_type" : "UserDefined",
"iter_max" : 50,
"delta" : 100,
"eps_trunc" : 1e-10,
"eps_dx" : 1e-16,
"glob_iter_max" : 30
}
}
This gives me:
ip_resolve('tr_newton.json')
iter f(x) ||grad|| ||dx|| mu_est
1 1.40e+01 6.39e+00 . 1.00e+00
2 4.12e+00 1.95e+00 1.72e+00 1.15e-01
. 4.12e+00 1.95e+00 7.29e-01 1.15e-01
3 5.58e+00 3.28e+01 3.64e-01 7.38e+00
4 5.32e+00 4.44e-16 1.61e-01 9.71e-01
5 3.95e+00 6.28e-16 3.20e-01 2.02e-01
6 3.76e+00 6.28e-16 5.29e-02 1.03e-01
7 3.58e+00 8.88e-16 4.67e-02 1.22e-02
8 3.56e+00 4.44e-16 5.94e-03 1.04e-03
9 3.56e+00 8.88e-16 4.96e-04 1.00e-04
10 3.56e+00 1.40e-15 4.79e-05 1.00e-05
11 3.56e+00 0.00e+00 4.77e-06 1.00e-06
12 3.56e+00 4.44e-16 4.77e-07 1.00e-07
13 3.56e+00 0.00e+00 4.77e-08 1.00e-08
The optimal point for c=1.000000e+00 is: (3.333333e-01,3.333333e-01)
iter f(x) ||grad|| ||dx|| mu_est
1 3.56e+00 0.00e+00 4.77e-08 1.00e-08
. 3.56e+00 0.00e+00 3.64e-01 1.00e-08
. 3.56e+00 0.00e+00 1.82e-01 1.00e-08
. 3.56e+00 0.00e+00 9.11e-02 1.00e-08
. 3.56e+00 0.00e+00 4.55e-02 1.00e-08
. 3.56e+00 0.00e+00 2.28e-02 1.00e-08
. 3.56e+00 0.00e+00 1.14e-02 1.00e-08
. 3.56e+00 0.00e+00 5.69e-03 1.00e-08
. 3.56e+00 0.00e+00 2.85e-03 1.00e-08
. 3.56e+00 0.00e+00 1.42e-03 1.00e-08
. 3.56e+00 0.00e+00 7.12e-04 1.00e-08
. 3.56e+00 0.00e+00 3.56e-04 1.00e-08
. 3.56e+00 0.00e+00 1.78e-04 1.00e-08
. 3.56e+00 0.00e+00 8.89e-05 1.00e-08
. 3.56e+00 0.00e+00 4.45e-05 1.00e-08
. 3.56e+00 0.00e+00 2.22e-05 1.00e-08
. 3.56e+00 0.00e+00 1.11e-05 1.00e-08
. 3.56e+00 0.00e+00 5.56e-06 1.00e-08
. 3.56e+00 0.00e+00 2.78e-06 1.00e-08
. 3.56e+00 0.00e+00 1.39e-06 1.00e-08
. 3.56e+00 0.00e+00 6.95e-07 1.00e-08
. 3.56e+00 0.00e+00 3.47e-07 1.00e-08
. 3.56e+00 0.00e+00 1.74e-07 1.00e-08
. 3.56e+00 0.00e+00 8.69e-08 1.00e-08
2 3.56e+00 3.77e+08 4.34e-08 9.88e+06
...
39 3.39e+00 9.93e-16 6.86e-05 1.00e-05
40 3.39e+00 4.44e-16 6.65e-06 1.00e-06
41 3.39e+00 4.44e-16 6.64e-07 1.00e-07
42 3.39e+00 6.28e-16 6.64e-08 1.00e-08
The optimal point for c=7.777778e-01 is: (4.074074e-01,1.851852e-01)
iter f(x) ||grad|| ||dx|| mu_est
1 3.39e+00 6.28e-16 6.64e-08 1.00e-08
. 3.39e+00 6.28e-16 1.82e-01 1.00e-08
. 3.39e+00 6.28e-16 9.11e-02 1.00e-08
. 3.39e+00 6.28e-16 4.55e-02 1.00e-08
. 3.39e+00 6.28e-16 2.28e-02 1.00e-08
. 3.39e+00 6.28e-16 1.14e-02 1.00e-08
. 3.39e+00 6.28e-16 5.69e-03 1.00e-08
. 3.39e+00 6.28e-16 2.85e-03 1.00e-08
. 3.39e+00 6.28e-16 1.42e-03 1.00e-08
. 3.39e+00 6.28e-16 7.12e-04 1.00e-08
. 3.39e+00 6.28e-16 3.56e-04 1.00e-08
. 3.39e+00 6.28e-16 1.78e-04 1.00e-08
. 3.39e+00 6.28e-16 8.89e-05 1.00e-08
. 3.39e+00 6.28e-16 4.45e-05 1.00e-08
. 3.39e+00 6.28e-16 2.22e-05 1.00e-08
. 3.39e+00 6.28e-16 1.11e-05 1.00e-08
. 3.39e+00 6.28e-16 5.56e-06 1.00e-08
. 3.39e+00 6.28e-16 2.78e-06 1.00e-08
. 3.39e+00 6.28e-16 1.39e-06 1.00e-08
. 3.39e+00 6.28e-16 6.95e-07 1.00e-08
. 3.39e+00 6.28e-16 3.47e-07 1.00e-08
. 3.39e+00 6.28e-16 1.74e-07 1.00e-08
. 3.39e+00 6.28e-16 8.69e-08 1.00e-08
2 3.39e+00 3.68e+08 4.34e-08 7.13e+06
...
37 3.48e+00 1.05e-14 1.03e-01 1.11e-01
38 3.30e+00 1.40e-15 6.37e-02 1.41e-02
39 3.27e+00 1.26e-15 2.56e-02 1.65e-03
40 3.27e+00 0.00e+00 4.01e-03 1.16e-04
41 3.27e+00 1.26e-15 1.97e-04 1.00e-05
42 3.27e+00 1.40e-15 1.33e-05 1.00e-06
43 3.27e+00 4.44e-16 1.32e-06 1.00e-07
44 3.27e+00 1.33e-15 1.32e-07 1.00e-08
The optimal point for c=5.555556e-01 is: (4.814815e-01,3.703705e-02)
iter f(x) ||grad|| ||dx|| mu_est
1 3.27e+00 1.33e-15 1.32e-07 1.00e-08
. 3.27e+00 1.33e-15 3.64e-01 1.00e-08
. 3.27e+00 1.33e-15 1.82e-01 1.00e-08
. 3.27e+00 1.33e-15 9.11e-02 1.00e-08
. 3.27e+00 1.33e-15 4.55e-02 1.00e-08
. 3.27e+00 1.33e-15 2.28e-02 1.00e-08
. 3.27e+00 1.33e-15 1.14e-02 1.00e-08
. 3.27e+00 1.33e-15 5.69e-03 1.00e-08
. 3.27e+00 1.33e-15 2.85e-03 1.00e-08
. 3.27e+00 1.33e-15 1.42e-03 1.00e-08
. 3.27e+00 1.33e-15 7.12e-04 1.00e-08
. 3.27e+00 1.33e-15 3.56e-04 1.00e-08
. 3.27e+00 1.33e-15 1.78e-04 1.00e-08
. 3.27e+00 1.33e-15 8.89e-05 1.00e-08
. 3.27e+00 1.33e-15 4.45e-05 1.00e-08
. 3.27e+00 1.33e-15 2.22e-05 1.00e-08
. 3.27e+00 1.33e-15 1.11e-05 1.00e-08
. 3.27e+00 1.33e-15 5.56e-06 1.00e-08
. 3.27e+00 1.33e-15 2.78e-06 1.00e-08
. 3.27e+00 1.33e-15 1.39e-06 1.00e-08
. 3.27e+00 1.33e-15 6.95e-07 1.00e-08
. 3.27e+00 1.33e-15 3.47e-07 1.00e-08
. 3.27e+00 1.33e-15 1.74e-07 1.00e-08
2 3.27e+00 3.62e+08 8.69e-08 4.39e+06
...
39 3.21e+00 5.07e-15 8.56e-04 1.07e-05
40 3.21e+00 6.04e-15 4.94e-05 1.00e-06
41 3.21e+00 8.01e-16 4.20e-06 1.00e-07
42 3.21e+00 2.09e-15 4.18e-07 1.00e-08
The optimal point for c=3.333333e-01 is: (5.555555e-01,-1.111111e-01)
iter f(x) ||grad|| ||dx|| mu_est
1 3.21e+00 2.09e-15 4.18e-07 1.00e-08
. 3.21e+00 2.09e-15 3.64e-01 1.00e-08
. 3.21e+00 2.09e-15 1.82e-01 1.00e-08
. 3.21e+00 2.09e-15 9.11e-02 1.00e-08
. 3.21e+00 2.09e-15 4.55e-02 1.00e-08
. 3.21e+00 2.09e-15 2.28e-02 1.00e-08
. 3.21e+00 2.09e-15 1.14e-02 1.00e-08
. 3.21e+00 2.09e-15 5.69e-03 1.00e-08
. 3.21e+00 2.09e-15 2.85e-03 1.00e-08
. 3.21e+00 2.09e-15 1.42e-03 1.00e-08
. 3.21e+00 2.09e-15 7.12e-04 1.00e-08
. 3.21e+00 2.09e-15 3.56e-04 1.00e-08
. 3.21e+00 2.09e-15 1.78e-04 1.00e-08
. 3.21e+00 2.09e-15 8.89e-05 1.00e-08
. 3.21e+00 2.09e-15 4.45e-05 1.00e-08
. 3.21e+00 2.09e-15 2.22e-05 1.00e-08
. 3.21e+00 2.09e-15 1.11e-05 1.00e-08
. 3.21e+00 2.09e-15 5.56e-06 1.00e-08
. 3.21e+00 2.09e-15 2.78e-06 1.00e-08
. 3.21e+00 2.09e-15 1.39e-06 1.00e-08
. 3.21e+00 2.09e-15 6.95e-07 1.00e-08
. 3.21e+00 2.09e-15 3.47e-07 1.00e-08
. 3.21e+00 2.09e-15 1.74e-07 1.00e-08
2 3.21e+00 3.58e+08 8.69e-08 1.65e+06
...
42 3.20e+00 1.60e-15 7.56e-05 1.01e-06
43 3.20e+00 1.25e-14 6.88e-06 1.00e-07
44 3.20e+00 7.77e-15 6.80e-07 1.00e-08
The optimal point for c=1.111111e-01 is: (6.000000e-01,-1.999999e-01)
iter f(x) ||grad|| ||dx|| mu_est
1 3.20e+00 7.77e-15 6.80e-07 1.00e-08
. 3.20e+00 7.77e-15 3.64e-01 1.00e-08
. 3.20e+00 7.77e-15 1.82e-01 1.00e-08
. 3.20e+00 7.77e-15 9.11e-02 1.00e-08
. 3.20e+00 7.77e-15 4.55e-02 1.00e-08
. 3.20e+00 7.77e-15 2.28e-02 1.00e-08
. 3.20e+00 7.77e-15 1.14e-02 1.00e-08
. 3.20e+00 7.77e-15 5.69e-03 1.00e-08
. 3.20e+00 7.77e-15 2.85e-03 1.00e-08
. 3.20e+00 7.77e-15 1.42e-03 1.00e-08
. 3.20e+00 7.77e-15 7.12e-04 1.00e-08
. 3.20e+00 7.77e-15 3.56e-04 1.00e-08
. 3.20e+00 7.77e-15 1.78e-04 1.00e-08
. 3.20e+00 7.77e-15 8.89e-05 1.00e-08
. 3.20e+00 7.77e-15 4.45e-05 1.00e-08
. 3.20e+00 7.77e-15 2.22e-05 1.00e-08
. 3.20e+00 7.77e-15 1.11e-05 1.00e-08
. 3.20e+00 7.77e-15 5.56e-06 1.00e-08
. 3.20e+00 7.77e-15 2.78e-06 1.00e-08
. 3.20e+00 7.77e-15 1.39e-06 1.00e-08
. 3.20e+00 7.77e-15 6.95e-07 1.00e-08
. 3.20e+00 7.77e-15 3.47e-07 1.00e-08
. 3.20e+00 7.77e-15 1.74e-07 1.00e-08
2 3.20e+00 3.58e+08 8.69e-08 1.78e+01
...
39 3.20e+00 9.16e-16 2.01e-05 1.00e-06
40 3.20e+00 6.28e-16 1.96e-06 1.00e-07
41 3.20e+00 2.09e-15 1.96e-07 1.00e-08
The optimal point for c=-1.111111e-01 is: (6.000000e-01,-2.000000e-01)
iter f(x) ||grad|| ||dx|| mu_est
1 3.20e+00 2.09e-15 1.96e-07 1.00e-08
. 3.20e+00 2.09e-15 3.64e-01 1.00e-08
. 3.20e+00 2.09e-15 1.82e-01 1.00e-08
. 3.20e+00 2.09e-15 9.11e-02 1.00e-08
. 3.20e+00 2.09e-15 4.55e-02 1.00e-08
. 3.20e+00 2.09e-15 2.28e-02 1.00e-08
. 3.20e+00 2.09e-15 1.14e-02 1.00e-08
. 3.20e+00 2.09e-15 5.69e-03 1.00e-08
. 3.20e+00 2.09e-15 2.85e-03 1.00e-08
. 3.20e+00 2.09e-15 1.42e-03 1.00e-08
. 3.20e+00 2.09e-15 7.12e-04 1.00e-08
. 3.20e+00 2.09e-15 3.56e-04 1.00e-08
. 3.20e+00 2.09e-15 1.78e-04 1.00e-08
. 3.20e+00 2.09e-15 8.89e-05 1.00e-08
. 3.20e+00 2.09e-15 4.45e-05 1.00e-08
. 3.20e+00 2.09e-15 2.22e-05 1.00e-08
. 3.20e+00 2.09e-15 1.11e-05 1.00e-08
. 3.20e+00 2.09e-15 5.56e-06 1.00e-08
. 3.20e+00 2.09e-15 2.78e-06 1.00e-08
. 3.20e+00 2.09e-15 1.39e-06 1.00e-08
. 3.20e+00 2.09e-15 6.95e-07 1.00e-08
. 3.20e+00 2.09e-15 3.47e-07 1.00e-08
. 3.20e+00 2.09e-15 1.74e-07 1.00e-08
2 3.20e+00 3.58e+08 8.69e-08 1.69e+01
...
39 3.20e+00 9.93e-16 1.16e-05 1.00e-06
40 3.20e+00 2.22e-16 1.16e-06 1.00e-07
41 3.20e+00 2.22e-16 1.16e-07 1.00e-08
The optimal point for c=-3.333333e-01 is: (6.000000e-01,-2.000000e-01)
iter f(x) ||grad|| ||dx|| mu_est
1 3.20e+00 2.22e-16 1.16e-07 1.00e-08
. 3.20e+00 2.22e-16 3.64e-01 1.00e-08
. 3.20e+00 2.22e-16 1.82e-01 1.00e-08
. 3.20e+00 2.22e-16 9.11e-02 1.00e-08
. 3.20e+00 2.22e-16 4.55e-02 1.00e-08
. 3.20e+00 2.22e-16 2.28e-02 1.00e-08
. 3.20e+00 2.22e-16 1.14e-02 1.00e-08
. 3.20e+00 2.22e-16 5.69e-03 1.00e-08
. 3.20e+00 2.22e-16 2.85e-03 1.00e-08
. 3.20e+00 2.22e-16 1.42e-03 1.00e-08
. 3.20e+00 2.22e-16 7.12e-04 1.00e-08
. 3.20e+00 2.22e-16 3.56e-04 1.00e-08
. 3.20e+00 2.22e-16 1.78e-04 1.00e-08
. 3.20e+00 2.22e-16 8.89e-05 1.00e-08
. 3.20e+00 2.22e-16 4.45e-05 1.00e-08
. 3.20e+00 2.22e-16 2.22e-05 1.00e-08
. 3.20e+00 2.22e-16 1.11e-05 1.00e-08
. 3.20e+00 2.22e-16 5.56e-06 1.00e-08
. 3.20e+00 2.22e-16 2.78e-06 1.00e-08
. 3.20e+00 2.22e-16 1.39e-06 1.00e-08
. 3.20e+00 2.22e-16 6.95e-07 1.00e-08
. 3.20e+00 2.22e-16 3.47e-07 1.00e-08
. 3.20e+00 2.22e-16 1.74e-07 1.00e-08
2 3.20e+00 3.58e+08 8.69e-08 1.67e+01
...
40 3.20e+00 6.28e-16 8.38e-07 1.00e-07
41 3.20e+00 4.97e-16 8.38e-08 1.00e-08
The optimal point for c=-5.555556e-01 is: (6.000000e-01,-2.000000e-01)
iter f(x) ||grad|| ||dx|| mu_est
1 3.20e+00 4.97e-16 8.38e-08 1.00e-08
. 3.20e+00 4.97e-16 3.64e-01 1.00e-08
. 3.20e+00 4.97e-16 1.82e-01 1.00e-08
. 3.20e+00 4.97e-16 9.11e-02 1.00e-08
. 3.20e+00 4.97e-16 4.55e-02 1.00e-08
. 3.20e+00 4.97e-16 2.28e-02 1.00e-08
. 3.20e+00 4.97e-16 1.14e-02 1.00e-08
. 3.20e+00 4.97e-16 5.69e-03 1.00e-08
. 3.20e+00 4.97e-16 2.85e-03 1.00e-08
. 3.20e+00 4.97e-16 1.42e-03 1.00e-08
. 3.20e+00 4.97e-16 7.12e-04 1.00e-08
. 3.20e+00 4.97e-16 3.56e-04 1.00e-08
. 3.20e+00 4.97e-16 1.78e-04 1.00e-08
. 3.20e+00 4.97e-16 8.89e-05 1.00e-08
. 3.20e+00 4.97e-16 4.45e-05 1.00e-08
. 3.20e+00 4.97e-16 2.22e-05 1.00e-08
. 3.20e+00 4.97e-16 1.11e-05 1.00e-08
. 3.20e+00 4.97e-16 5.56e-06 1.00e-08
. 3.20e+00 4.97e-16 2.78e-06 1.00e-08
. 3.20e+00 4.97e-16 1.39e-06 1.00e-08
. 3.20e+00 4.97e-16 6.95e-07 1.00e-08
. 3.20e+00 4.97e-16 3.47e-07 1.00e-08
. 3.20e+00 4.97e-16 1.74e-07 1.00e-08
2 3.20e+00 3.58e+08 8.69e-08 1.67e+01
...
39 3.20e+00 0.00e+00 6.75e-05 1.00e-05
40 3.20e+00 9.93e-16 6.67e-06 1.00e-06
41 3.20e+00 4.44e-16 6.67e-07 1.00e-07
42 3.20e+00 4.97e-16 6.67e-08 1.00e-08
The optimal point for c=-7.777778e-01 is: (6.000000e-01,-2.000000e-01)
iter f(x) ||grad|| ||dx|| mu_est
1 3.20e+00 4.97e-16 6.67e-08 1.00e-08
. 3.20e+00 4.97e-16 8.55e-02 1.00e-08
. 3.20e+00 4.97e-16 4.28e-02 1.00e-08
. 3.20e+00 4.97e-16 2.14e-02 1.00e-08
. 3.20e+00 4.97e-16 1.07e-02 1.00e-08
. 3.20e+00 4.97e-16 5.35e-03 1.00e-08
. 3.20e+00 4.97e-16 2.67e-03 1.00e-08
. 3.20e+00 4.97e-16 1.34e-03 1.00e-08
. 3.20e+00 4.97e-16 6.68e-04 1.00e-08
. 3.20e+00 4.97e-16 3.34e-04 1.00e-08
. 3.20e+00 4.97e-16 1.67e-04 1.00e-08
. 3.20e+00 4.97e-16 8.35e-05 1.00e-08
. 3.20e+00 4.97e-16 4.18e-05 1.00e-08
. 3.20e+00 4.97e-16 2.09e-05 1.00e-08
. 3.20e+00 4.97e-16 1.04e-05 1.00e-08
. 3.20e+00 4.97e-16 5.22e-06 1.00e-08
. 3.20e+00 4.97e-16 2.61e-06 1.00e-08
. 3.20e+00 4.97e-16 1.31e-06 1.00e-08
. 3.20e+00 4.97e-16 6.53e-07 1.00e-08
. 3.20e+00 4.97e-16 3.26e-07 1.00e-08
. 3.20e+00 4.97e-16 1.63e-07 1.00e-08
2 3.20e+00 3.58e+08 8.16e-08 1.57e+01
3 3.20e+00 3.58e+06 8.16e-08 6.57e-01
...
38 3.20e+00 2.22e-16 6.01e-04 1.00e-04
39 3.20e+00 4.97e-16 5.67e-05 1.00e-05
40 3.20e+00 0.00e+00 5.63e-06 1.00e-06
41 3.20e+00 4.97e-16 5.63e-07 1.00e-07
42 3.20e+00 4.97e-16 5.63e-08 1.00e-08
The optimal point for c=-1.000000e+00 is: (6.000000e-01,-2.000000e-01)
Note, I truncated things since I ran over the character limit for this board. Alright, so it works, but not exceptionally well. Let me explain why. Basically, the inequalities are handled with a primal dual interior point method. Generally, the name of the game in such methods is to play games with the barrier parameter mu to help guide the solution into the optimal point. If we run into the boundary too quickly, the algorithm may be able to recover, but generally it’s going to be slow since it’s going to run itself into the boundary over and over again. There is a body of research that talks about how to warm start an interior point method effectively, but generally that work has been done with linear, rather than nonlinear, cone problems. Perhaps it carries over well, but I’ve never investigated it closely.
Alright, none of that really helps you all that much. Here’s the immediate fix:
- Set
opt_stop to NotConverged
- Set
iter to 1
- Set
mu to something kind of big to help push off the boundary. In the absence of a good idea, I use 1.
- Set
glob_iter_max to something big because the first iteration is going to be rough.
Now, are there better options than this? Yes. Personally, if you’re going to do a lot of warm starts and the solution doesn’t change that much, I like the reflective Newton method from Coleman and Li in their paper “An Interior Trust Region Approach For Nonlinear Minimization Subject to Bounds”. I find the reflective stuff to be kind of a gimmick. That said, if your solutions are close together so that you stay inside the Newton convergence radius, then this algorithm is really good and it doesn’t care if the new problem makes the old problem infeasible. If the solutions are farther away, then things may get more dicey.
If your problem was linear, then the easy answer would be to just add the new constraint and solve the dual problem with the simplex method (or, really, dual simplex.) If the solution is close by, then perhaps only a couple of pivots are necessary, which are fast. This is what integer programming solvers do when implementing branch and bound. For nonlinear problems, we most often can’t play the duality trick. Maybe there’s a trick. Likely an active set method would work fine, but there probably needs to be some care in determining the initial set of active constraints and dealing with adding a new constraint that breaks feasibility.
Alright, so there are other algorithms that work for this case better. That said, unfortunately, Optizelle doesn’t implement them! That’s, well, unfortunate. I am working on an updated version of the interior point method that would perhaps make this easier, but, really, a primal-dual interior point method isn’t super well suited for this kind of task unless there’s been some new update in research that I’m not aware of.
Now, at the very least, the code should not crash. That’s worrying. If you have something manageably sized that demonstrates the crash, I’d like to have a look.
Short of that, use the trick above to maybe, sort of, force it to work. Ultimately though, you’ll likely need a different algorithm.
