fmincon stopped because the size of the current step is less than...

I apologize for the lengthly code below, but if someone can figure out the problem I will be gateful!

I receive th efollowing message for this code:

Converged to an infeasible point.

fmincon stopped because the size of the current step is less than

the default value of the step size tolerance but constraints are not

satisfied to within the selected value of the constraint tolerance.

<stopping criteria details>

clear all;
mpc = case14;
%% data for the system
p=zeros(14,1);
p(1,1)=220.97;
p(2,1)=38.03;
d=mpc.bus(:,3); % Actual value of load d(in MW)
G=makePTDF(case14);%Shifting factor matrix
% left null space
[V, pivot] = rref(G)
r = length(pivot)
cs = G(:,pivot)
ns = null(G,'r')
rs = V(1:r,:)'
lns = null(G','r')
H=lns'
B=pinv(G);
c=repmat(39.016,14,1);
delta_p_low=zeros(14,1);
for i=1:14
    delta_p_low(i,1)=-2;
end
delta_p_up=zeros(14,1);
for i=1:6
    delta_p_up(i,1)=0.1;
end
f_up=repmat(200,20,1);
f_low=repmat(-200,20,1);
%% variables
 lamda_plus=sdpvar(14,1);
 lamda_minus=sdpvar(14,1);
 alpha=sdpvar(14,1);
 f=G*(p-alpha-d);
 miu_minus=sdpvar(20,1);  
 miu_plus=sdpvar(20,1);
 v=sdpvar(7,1); 

%% constraints
F=[];
F=[F 0<=alpha(1)<=p(1)];
F=[F alpha(2,1)==0];
F=[F alpha(3,1)==0];
F=[F alpha(4,1)==0];
F=[F alpha(5,1)==0];
F=[F alpha(6,1)==0];
F=[F alpha(7,1)==0];
F=[F alpha(8,1)==0];
F=[F alpha(9,1)==0];
F=[F alpha(10,1)==0];
F=[F alpha(11,1)==0];
F=[F alpha(12,1)==0];
F=[F alpha(13,1)==0];
F=[F alpha(14,1)==0];
F=[F delta_p_low<=B*f-p+alpha+d<=delta_p_up];
F=[F f_low<=f<=f_up];
F=[F H*f==0];
F=[F lamda_minus>=0];
F=[F lamda_plus>=0];
F=[F miu_minus>=0];
F=[F miu_plus>=0];
F=[F diag(lamda_plus*((B*f)-p+alpha+d-delta_p_up)')==zeros(14,1)];
F=[F diag(lamda_minus*(delta_p_low-(B*f)+p-alpha-d)')==zeros(14,1)];
F=[F diag(miu_plus*(f-f_up)')==zeros(20,1)];
F=[F diag(miu_minus*(f_low-f)')==zeros(20,1)];
F=[F B'*(c+lamda_plus-lamda_minus)+miu_plus-miu_minus+H'*v==0];
%% objective
obj=-((c(1)+lamda_plus(1)-lamda_minus(1))*(p(1)-alpha(1))-c(1)*p(1));
%% solution
% option=optimset( 'Algorithm', 'interior-point','MaxFunEvals', 200000,'MaxIter',5000,'TolFun',1e-6,'TolX',1e-12);
option=sdpsettings('savesolveroutput',1);
diagnostics=optimize(F,obj,option);
result.obj=double(obj);

[Walter Roberson]

Is case14 https://www.mathworks.com/matlabcentral/fileexchange/43139-14bus

[Vahap Demir]

Thank you Walter, I appreciated. But I have already put the data into my code. I wonder what the mistake is in my code.

[Walter Roberson]

We do not have case14 . Is that particular place where your case14 is from? We cannot run your code without the data.

[Vahap Demir]

I used Case14 in matpower. But even with the data the code does not work. I appreciate any further help!

[Walter Roberson]

(Your code also requires YALMIP)

[Vahap Demir]

I am using YALMIP too. The problem persists! Is there anything else you can suggest?

Thank you again!

[Walter Roberson]

Note: to get the condition you mention, I had to change the last bit of your code to

%% solution
fmincon_option=optimset( 'Algorithm', 'interior-point','MaxFunEvals', 200000,'MaxIter',50000,'TolFun',1e-6,'TolX',1e-12);
option=sdpsettings('savesolveroutput',1);
option.fmincon = fmincon_option;
diagnostics=optimize(F,obj,option)
result.obj=double(obj);

[Vahap Demir]

I really appreciated your time! After using your solution the problem still persists. I still receive the same error message. Please see below:

Answer 1

[Walter Roberson]

With the combinations of configuration and options you have used, yalmip ends up using all zeros as the initial point for fmincon. When starting from that point, fmincon ends up in a local minima in which not all of the constraints can be met.

Your situation has two possibilities:

1. Perhaps it is not possible to meet your constraints.
2. Perhaps with those fmincon settings and those initial conditions, it is not possible to reach a location that meets the constraints.

You have nonlinear constraints, including nonlinear equality constraints. A violation of those constraints, especially the equality ones, can be difficult to remedy, as the path to reach a valid point might be a curved line in 89 dimensional space. When the constraints can be met, a fortunate starting point can be crucial.

yalmip does not give a direct method of choosing a starting point. However, for each set of variables you define with sdpvar, you can use assign() https://yalmip.github.io/command/assign/ to give an initial value. Then in your sdpsettings call, also add 'usex0', 1 to the options. yalmip will use the values assigned to construct the x0 vector.

You might need to try a few different possibilities... many of them even. Finding a feasible point can be hard when it exists at all.

Answer 2

[Johan Löfberg]

Note that your model contains both a lot of numerical noise (things that you probably think is zero, but is on the order of 10^-15 and smaller), and a completely wird 0*f = 0 constraint

The 17th constraint corresponds to H*f == 0, which ends up being almost 0 == 0, but not really due to some remaining numerical noise

F([15 17 22 26])
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|                            Constraint|          Coefficient range|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|          Element-wise inequality 28x1|       3.4694e-18 to 221.07|
|   #2|               Equality constraint 7x1|   8.6736e-19 to 7.1054e-14|
|   #3|   Equality constraint (bilinear) 14x1|       1.7764e-14 to 221.07|
|   #4|              Equality constraint 20x1|      2.7933e-19 to 44.7174|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Hence, I would clean up B and nls before using them, and don't add H*f == 0 since it already holds in your model as H*f is 0 by construction

lns(abs(lns)<1e-12)=0;
B(abs(B)<1e-14)=0;

BTW, your code is unecessarily long, you can simply write

F = [F, alpha(2:end)==0]

Of course, one wonders why you even define all those variables when they are 0, instead of simply defining

alpha = [sdpvar(1);zeros(13,1)]

And why constructing a large symbolic square matrix

lamda_plus*((B*f)-p+alpha+d-delta_p_up)'

and then taking the diagonal of that. You simply want

lamda_plus.*((B*f)-p+alpha+d-delta_p_up)

same on all those final expressions you have in your complementary slackness conditions

Answer 3

[Vahap Demir]

Thank you very much! I really appreciated!