fmincon: any way to enforce linear inequality constraints at intermediate iterations?

I am solving an optimization problem with the interior-point algorithm of fmincon.

My parameters have both lower and upper bounds and linear inequality constraints.

A_ineq = [1 -1 0 0 0;
          -1 2 -1 0 0;
          0 -1 2 -1 0;
          0 0 -1 2 -1;
          0 0 0 1 -1];
b_ineq=[0 0 0 0 0];

I observed that fmincon satisfies, of corrse, the bounds at all iterations but not the linear inequality constraints.

However, the linear inequality constraints are crucially important in my case. If violated, my optimization problem can not be solved.

Is there anything one can do to ensure that linear inequality constraints are satisfied at intermediate iterations?

[Torsten]

No.

[SA-W]

Here, Alan Weiss pointed out to multiple the constraint function with an large value. Is it reasonable to multiply A_ineq by, say, 1e4 to ensure a better consideration of the inequality constraints?

[Walter Roberson]

https://www.mathworks.com/help/optim/ug/iterations-can-violate-constraints.html

fmincon interior-point, sqp, trust-region-reflective satisfy bounds constraints at each iteration, but active-set can violate bounds constraints.

But those refer to bound constraints, not to linear constraints.

[SA-W]

@Walter Roberson

I am aware of that. But I thought that there are maybe some manipulationens to the linear constraints such that fmincon gives them a higher prority at intermediate iterations.

[SA-W]

Scaling the matrix A_ineq by 1e4 or any other large value really helps in my case to (in a least-squares sense) enforce the linear constraints at intermediate iterations.

[Matt J]

I'd like to point out that it violates the theoretical assumptions of fmincon, and probably all the Optimization Toolbox solvers as well, when the domain of your objective function and constraints is not an open subset of . If you forbid evaluation outside the closed set defined by your inequality constraints, that is essentially the situation you are creating. It's not entirely clear what hazards this creates in practice, but it might mean one of the Global Optimization Toolbox solvers might be more appropriate.

Answer 1

[Matt J]

Within your objective function, project the current x onto the constrained region using quadprog,

fun=@(x) myObjective(x, A_ineq,b_ineq,lb,ub);

x=fmincon(fun, x0,[],[],[],[],[],lb,ub,options)

function fval=myObjective(x, A_ineq,b_ineq,lb,ub)

  C=speye(numel(x));
  x=lsqlin(C,x,A_ineq,b_ineq,[],[],lb,ub); %EDIT: was quadprog by mistake
  
  fval=...

end

[SA-W]

If I understand this snippet correctly, you want to invoke fmincon without passing any constraints and bounds and, at every iteration of fmincon, let quadprog project the x such that all constraints are satisfied?

If so I have two questions:

If constraints were satisfied, quadprog would not change the x, right?

I could also use this strategy with lsqnonlin, right?

[Matt J]

I changed my mind about passing bounds. You should put those in. Yes, I guess you could do it with lsqnonlin as well.

I'm not sure about the performance of this strategy, as compared to my other answer. It violates differentiability assumptions. Also, because you will be using finite difference derivative approximation, quadprog will need to be called N times per iteration, so it wouldn't be a good strategy for high-dimensional problems.

[SA-W]

And why would you not pass the linear constraints to fmincon but only to quadprog?

[Matt J]

Because the projection problem is quadratic.

[SA-W]

Can you elaborate a little bit more on that?

[Matt J]

It is a quadratic minimization problem

where C is the set defined by the constraints.

[SA-W]

I know. But why should I not pass the linear constraints to fmincon? I do not see the connection with the subsequent call to quadprog.

[Matt J]

Ah, well there's no reason fmincon should need to worry about enforcing linear inequality constraints if every test point is going to be projected into a feasible candidate anyway. The objective function is incapable of evaluating anything outside C.

Come to think of it, though, it would probably be wise to apply a penalty on distance from the feasible set as well to remove ambiguity in the solution.

function fval=myObjective(x, A_ineq,b_ineq,lb,ub, weight )

  C=speye(numel(x));
  [x,dist]=quadprog(C,x,A_ineq,b_ineq,[],[],lb,ub);
  
  fval=...
      
  fval=fval+weight*dist
      
  
end

[SA-W]

I tried your suggested strategy, but it is not working as expected:

function fval=myObjective(x, A_ineq,b_ineq,lb,ub, weight )
  
  %1st iteration, i.e. x is the start vector
  
  x = 
       0.004000000000000
       0.001000000000000
                       0
       0.001000000000000
       0.004000000000000
       0.009000000000000
       0.016000000000000
       0.025000000000000
       0.036000000000000
                       0
       0.020250000000000
       0.063245553203368
       0.083495553203368
       0.081000000000000
       0.144245553203368
       0.089442719099992
       0.109692719099992
       0.170442719099992
       
  all(A_ineq*x) % == 1, i.e. start vector satisfies constraints
  
  C=speye(numel(x));
  [x,dist]=quadprog(C,x,A_ineq,b_ineq,[],[],lb,ub);
  
  x =
      1.0e-03 *
    
       0.178567411241192
       0.050029327861103
                       0
       0.000014275626508
       0.000030351196081
       0.000053745685964
       0.000087040217878
       0.000137579529376
       0.000245546982784
                       0
       0.000003261476952
       0.000013844333853
       0.000023204570239
       0.000023729908401
       0.000054189214599
       0.000019086149826
       0.000029710389433
       0.000063797849365
       
   % not possible to work with the new x, which is several magnitudes lower
      
  
end

Calling quadprog with the start vector, which satisfies all constraints and bounds, results in a new vector x, where nearly all components are several orders of magnitudes lower.

I attached the matrix A_ineq, b_ineq=0, lb=2, ub=0. Basically, A_ineq encodes that parameters are less or greater than neighboring parameters.

Ideally, if constraints were satisfied, I would expect quadprog to not change x significantly.

Does it make sense to you how the projected x looks like in my case?

[SA-W]

Update:

The issue seems to be not related to the linear constraints but to the bounds.

C=speye(numel(x));
  [x,dist]=quadprog(C,x,[],[],[],[],lb,ub);

x = 
   1.0e-03 *

   0.001101980588430
   0.160659227282455
                   0
   0.160659227282455
   0.001101980588430
   0.000065909263473
   0.000037138227992
   0.000023769741153
   0.000016507062926
                   0
   0.000029344843098
   0.000009396094623
   0.000007117296309
   0.000007336573768
   0.000004119807053
   0.000006644060468
   0.000005417528065
   0.000003486590875
  
lb = 
     0
     0
     0
     0
     0
     0
     0
     0
     0
     0
     0
     0
     0
     0
     0
     0
     0
     0
   
ub = 
     1
     1
     0
     1
     1
     1
     1
     1
     1
     0
     2
     2
     2
     2
     2
     2
     2
     2

Calling quadprog without bounds and linear constraints, I get the the vector that I pass multiplied by (-1).

@Matt J

Do you have any idea as to why the bounds are causing this issue?

[Torsten]

You want to find a vector x that minimizes

1/2*(x-x_passed).'*(x-x_passed)

under the constraints

A_ineq*x <= b_ineq

where x_passed is the vector you get from fmincon.

Thus f in the list of inputs to quadprog must be -x_passed, not x_passed:

[x,dist]=quadprog(C,-x,A_ineq,b_ineq,[],[],lb,ub);

instead of

[x,dist]=quadprog(C,x,A_ineq,b_ineq,[],[],lb,ub);

[Matt J]

@SA-W Sorry, I made a mistake. Instead of quadprog, I meant for you to use lsqlin:

fun=@(x) myObjective(x, A_ineq,b_ineq,lb,ub);

x=fmincon(fun, x0,[],[],[],[],[],lb,ub,options)

function fval=myObjective(x, A_ineq,b_ineq,lb,ub)

  C=speye(numel(x));
  x=lsqlin(C,x,A_ineq,b_ineq,[],[],lb,ub);
  
  fval=...

end

But I want to reiterate that I think this whole thing you are pursuing is probably very inefficient. You should just abort the objective function with inf, as I suggest in my other answer.

[SA-W]

Thanks to both of you. Makes sense to pass -x to quadprog.

@Matt J I have an analytical expression for the Jacobian. Hence, quadprog of lsqlin have to be called only once per iteration of lsqnonlin/fmincon. That said, I think the overall strategy is not so inefficient.

Anyway, why would you use lsqlin instead of quadprog?

[Matt J]

@SA-W In the code I originally posted, quadprog was being called with the input argument syntax of lsqlin, so that probably explains why it didn't work. lsqlin is more specialized than quadprog, so I tend to think it is better to use that when it is applicable.

Regarding your analytical gradient, I am not sure it is valid anymore. If your original objective function was f(x) with gradient g(x), the new one with the projection onto the linear contraints is f(P(x)) where P is my notation for the projection operator. The new gradient would be *g(P(x)), but I don't immediately see any easy way to compute .

EDIT: Possibly, you can take as the linear projection operator (transposed) onto the active linear constraints.

[SA-W]

"The new gradient would be *g(P(x)), but I don't immediately see any easy way to compute nabla P "

I see what you mean. If I were using finite differencing at fmincon/lsqnonlin level, I would not have to worry about that, right?

" EDIT: Possibly, you can take nabla P as the linear projection operator (transposed) onto the active linear constraints. "

Does this help in calculating nabla P?

[Matt J]

I see what you mean. If I were using finite differencing at fmincon/lsqnonlin level, I would not have to worry about that, right?

Right, but it would slow you down a lot.

Does this help in calculating nabla P?

Yes, projecting onto the active constraints is a simple matrix multiplication.

[SA-W]

"Yes, projecting onto the active constraints is a simple matrix multiplication."

Could you please explain a little bit how that could be implemented?

[Matt J]

The orthogonal projection of z onto the set {x|A*x=b} is

P(z)=z-pinv(A)*(A*z-b)

See also,

https://math.stackexchange.com/questions/1320363/orthogonal-projection-of-z-onto-the-affine-set-left-x-mid-a-x-b-rig

So, the gradient matrix of P would be,

nablaP =(I-pinv(A)*A).';

[SA-W]

Do the expressions look the same for inequailities instead of equalities, i.e. if the set is given by

set {x|A*x <= b} instead of set {x|A*x=b} ?

[Matt J]

No, but the case {x|A*x <= b} doesn't apply here. You will be projecting onto the active constraints. Active constraints are the ones satisfied with equality at the solution lsqlin gives you..

[SA-W]

Say A*x > 0 at some components with x given by fmincon. After invoking lsqlin, I observed that A*x=0 at those components. Is this the expected behavior?

[Matt J]

I don't know what you mean by "with x given by fmincon". In the approach we are talking about, fmincon can does not act independently of lsqlin. You are calling lsqlin inside the objective function, which is in turn called iteratively by fmincon.

[SA-W]

"I don't know what you mean by "with x given by fmincon". In the approach we are talking about, fmincon can does not act independently of lsqlin. You are calling lsqlin inside the objective function, which is in turn called iteratively by fmincon."

Yes, that's clear. Within the objective function, say some constraints are violated, i.e. A_ineq*x > 0 for some components. Then I call lsqlin to project x onto the feasible set and what I observe then is that A_ineq*x = 0 at those components, i.e. those constraints are active now.

Does it make sense that the constraints are active after calling lsqlin?

[Torsten]

They can be active or inactive after calling "lsqlin". The purpose of the call to "lsqlin" is to make the x from "fmincon" feasible, and feasible means: a_i*x = b_i or a_i*x < b_i.

[SA-W]

Suppose a_i*x > b_i before calling lsqlin. After calling lsqlin, is this constraint active or inactive, or is there no general statement?

[Matt J]

The general statement is this - an inequality constraint a_i*x <= b_i can either be active, inactive, or violated at a given point x. In ideal mathematical terms, these cases are as follows:

a_i*x < b_i (inactive)

a_i*x = b_i (active)

a_i*x > b_i (violated)

An equality constraint will always be active, when satisfied. Otherwise, it is violated.

The optimization strategy we are pursuing here requires that you find which constraints are active at the projected point generated by lsqlin. As you do this, you have to be mindful that you are working with a finite precision computer, and so you have to modify the ideal notion of what it means for a constraint to be active. You must instead consider a constraint active if,

abs( a_i*x - b_i ) < tol

where tol>=0 is a tolerance parameter to be chosen by you.

[SA-W]

The new gradient would be *g(P(x)), but I don't immediately see any easy way to compute nabla P .

EDIT: Possibly, you can take nabla P as the linear projection operator (transposed) onto the active linear constraints.

So, in the equation

nablaP =(I-pinv(A)*A).';

the matrix A contains only rows associated with active constraints?

[Matt J]

That's right.

[Torsten]

Suppose a_i*x > b_i before calling lsqlin. After calling lsqlin, is this constraint active or inactive, or is there no general statement?

No general statement is possible in my opinion.

Answer 2

[Shubham]

Hi SA-W,

Yes, there are a few options you can try to ensure that the linear inequality constraints are satisfied at intermediate iterations of the interior-point algorithm in fmincon:

1. Tighten the tolerances: By tightening the tolerances in the fmincon options, you can force the algorithm to take smaller steps and converge more slowly, but with higher accuracy. This may help to ensure that the linear inequality constraints are satisfied at all intermediate iterations.
2. Use a barrier function: You can try using a barrier function to penalize violations of the linear inequality constraints. This can be done by adding a term to the objective function that grows very large as the constraints are violated. This will encourage the algorithm to stay within the feasible region defined by the constraints.
3. Use a penalty function: Similar to a barrier function, a penalty function can be used to penalize violations of the linear inequality constraints. However, instead of growing very large, the penalty function grows linearly with the degree of violation. This can be a more computationally efficient approach than a barrier function.
4. Use a combination of methods: You can try using a combination of the above methods to ensure that the linear inequality constraints are satisfied at all intermediate iterations. For example, you could tighten the tolerances and use a penalty function or barrier function to further enforce the constraints.

It's important to note that these methods may increase the computational cost of the optimization problem, so it's important to balance the accuracy requirements with the available resources.

Answer 3

[Matt J]

Within your objective function, check if the linear inequalites are satsfied. If they are not, abort the function and return Inf. Otherwise, compute the output value as usual.

fun=@(x) myObjective(x, A_ineq,b_ineq);

x=fmincon(fun, x0,A_ineq,b_ineq,[],[],lb,ub,[],options)

function fval=myObjective(x, A_ineq,b_ineq)

  if ~all(A_ineq*x<=b_ineq)
     fval=Inf; return
  end
      
  fval=...

end