Is fmincon correct for this problem?

Hi, I am doing a power flow optimization. I am unable to solve it in MATLAB and perhaps the solver I am using (fmincon) is not the correct one. Guidance will be highly appreciated.

The problem is as follows:

% Defining my lower and upper bounds
lb_Vm_ac =  0.95*ones(4,1);
lb_Va_ac = -pi*ones(4,1);
lb_Pdg_ac = [-Inf; 0.1; 0.05];
lb_Qdg_ac = [-Inf; -0.2; -0.2];


ub_Vm_ac = 1.05*ones(4,1);
ub_Va_ac = pi*ones(4,1);
ub_Pdg_ac = [Inf; 0.4; 0.4];
ub_Qdg_ac = [Inf; 0.3; 0.2];


A = []; b = []; 
Aeq = []; beq = [];
lb = [lb_Vm_ac; lb_Va_ac; lb_Pdg_ac; lb_Qdg_ac]; 
ub = [ub_Vm_ac; ub_Va_ac; ub_Pdg_ac; ub_Qdg_ac];

% Defining Initial Guess

Vm_ac = 1*ones(4,1);
Va_ac = zeros(4,1);
Pdg_ac = 0.2*ones(3,1); 
Qdg_ac = 0.2*ones(3,1);

X0 = [Vm_ac; Va_ac; Pdg_ac; Qdg_ac]

% Executing Solution

[x,feval] = fmincon(@objective_loss,X0,A,b,Aeq,beq,lb,ub,@nl_constraints)

% Objective Function

function cost = objective_loss(X)
    Pdg_ac = X(9:11,1);

cost = 0.35*Pdg_ac(1,1) + 0.2*Pdg_ac(2,1) + 0.4*Pdg_ac(2,1)^2 + 0.3*Pdg_ac(3,1) + 0.5*Pdg_ac(3,1)^2;

% Nonlinear Constraints

function [c,ceq] = nl_constraints(X)
    Vm_ac = [1; X(1:4,1)];
    Va_ac = [0; X(5:8,1)];
    Pdg_ac = X(9:11,1);
    Qdg_ac = X(12:14,1);
    
    Pg_ac = [Pdg_ac(1,1); 0; Pdg_ac(2,1); Pdg_ac(3,1); 0];
    Pd_ac = [0; 0; 0; 0.9; 0.239];
    Qg_ac = [Qdg_ac(1,1); 0; Qdg_ac(2,1); Qdg_ac(3,1); 0];
    Qd_ac = [0; 0; 0; 0.4; 0.129];
    
    Pac = (Pg_ac - Pd_ac); 
    Qac = (Qg_ac - Qd_ac);
    
    
       
    Y = [1.07-10.04i 0+3.33i -1.07+6.73i 0 0;
        0+3.33i 5.66-33.33i 0 -5.66+30.19i 0;
        -1.07+6.73i 0 1.41-13.78i -0.09+3.83i 0+3.19i;
        0 -5.66+30.19i -0.49+3.8i 5.95-36.01i 0+2i;
        0 0 0+3.19i 0+2i 0-5.13i];
    
    % Inequality Constraints
    Pgen_more_than_Pload = sum(Pd_ac) - sum(Pg_ac);
    Qgen_more_than_Qload = sum(Qd_ac) - sum(Qg_ac);
    
    
    
    % Equality Constraints
    power_flow1 = Pac - real(diag(Vm_ac.*exp(1i*Va_ac))*conj(Y*Vm_ac.*exp(1i*Va_ac)));
    power_flow2 = Qac - imag(diag(Vm_ac.*exp(1i*Va_ac))*conj(Y*Vm_ac.*exp(1i*Va_ac)));
    
c = [Pgen_more_than_Pload; Qgen_more_than_Qload];
ceq = [power_flow1; power_flow2];

I am not getting a feasible solution. This is what I get

>> [x,feval] = fmincon(@objective_loss,X0,A,b,Aeq,beq,lb,ub,@nl_constraints)

Converged to an infeasible point.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance but constraints are not
satisfied to within the value of the constraint tolerance.

<stopping criteria details>

x =

    0.9588
    0.9923
    0.9500
    0.9577
   -3.1346
   -3.1413
    0.0691
    0.9723
    0.2217
    0.4000
    0.4000
    0.2382
    0.3000
   -0.0476


feval =

    0.4216

However, I am sure this is solvable. As I am just coding this problem from a manuscript.

Kindly advise. Please note that some of my constraints use real(), imag() finctions. As this is a power flow problem, it has complex quantities.

[Matt J]

However, I am sure this is solvable. As I am just coding this problem from a manuscript.

Is the choice of X0 from the manuscript as well?

[Arif Ahmed]

The manuscript doesn't mention the choice of X0 as they used GAMS, which I believe can solve without an initial estimate.

[Matt J]

which I believe can solve without an initial estimate.

No, that would not be possible. You should contact the authors and ask their advice on initialization.

Answer 1

[Alan Weiss]

I see that you define a quantity Y in your nonlinear constraint expression, but you don't seem to use it. Was this an oversight?

There are suggestions in the documentation for steps to take when you converge to an infeasible point.

Alan Weiss

MATLAB mathematical toolbox documentation

[Arif Ahmed]

Thanks Alan,

Y is being used in power_flow1 and power_flow2