You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
fmincon: any way to enforce linear inequality constraints at intermediate iterations?
13 views (last 30 days)
Show older comments
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?
6 Comments
Walter Roberson
on 2 Mar 2023
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
on 2 Mar 2023
@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
on 2 Mar 2023
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
on 6 May 2023
Edited: Matt J
on 6 May 2023
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1376369/image.png)
Accepted Answer
Matt J
on 5 May 2023
Edited: Matt J
on 23 May 2023
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
29 Comments
SA-W
on 6 May 2023
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
on 6 May 2023
Edited: Matt J
on 6 May 2023
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
on 7 May 2023
And why would you not pass the linear constraints to fmincon but only to quadprog?
Matt J
on 7 May 2023
It is a quadratic minimization problem
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1376669/image.png)
where C is the set defined by the constraints.
SA-W
on 7 May 2023
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
on 7 May 2023
Edited: Matt J
on 7 May 2023
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
on 23 May 2023
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
on 23 May 2023
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).
Do you have any idea as to why the bounds are causing this issue?
Torsten
on 23 May 2023
Edited: Torsten
on 23 May 2023
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
on 23 May 2023
@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
on 23 May 2023
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
on 23 May 2023
Edited: Matt J
on 24 May 2023
@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
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1391479/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1391479/image.png)
EDIT: Possibly, you can take
as the linear projection operator (transposed) onto the active linear constraints.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1391479/image.png)
SA-W
on 23 May 2023
"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
on 23 May 2023
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
on 23 May 2023
"Yes, projecting onto the active constraints is a simple matrix multiplication."
Could you please explain a little bit how that could be implemented?
SA-W
on 23 May 2023
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
on 23 May 2023
Edited: Matt J
on 23 May 2023
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
on 23 May 2023
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?
SA-W
on 23 May 2023
"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
on 23 May 2023
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
on 24 May 2023
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
on 24 May 2023
Edited: Matt J
on 24 May 2023
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
on 24 May 2023
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?
Torsten
on 24 May 2023
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.
More Answers (2)
Shubham
on 5 May 2023
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:
- 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.
- 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.
- 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.
- 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.
Matt J
on 5 May 2023
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
See Also
Categories
Find more on Linear Least Squares in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)