I'm having trouble solving a system of nonlinear equations with the fmincon function.

I want to solve a system of equations containing two non-linear equations, one of which is Psh_newhv == Pse_newhv,the other is Qs_newhv==0, but solving it with the fmincon function doesn't achieve the first constraint, how can I solve this problem?

% Given constants
V_s = 10*1e3/sqrt(3);
n_t = 25*sqrt(3);
V_r = 400/sqrt(3);
theta_r = pi/12;
Srate = 100*1e3/3;
Zbase = (V_r)^2/Srate;
X_r = 0.4*Zbase;  
X_t = 400^2/(300*1e3)*4/100; 
rho = 0:2*pi/72:2*pi; 
V_se = 0.2*V_r; 
V_se_newhv = 0.2*V_s; 
Z_r = 1i * X_r;
% Initial guess
x0 = [0 0];
% Options for fmincon
options = optimoptions('fmincon', 'Algorithm', 'interior-point', 'Display', 'off','ConstraintTolerance',1e-6);
% Constraints
lb = [];
ub = [];
% Solve for each rho
for m = 1:73
    % Define the nonlinear constraint function
    fun_newhv = @(x_newhv) myfun_newhv(x_newhv, rho(m));
    % Use fmincon to solve the problem
    [x_newhv, fval_newhv, exitflag_newhv(m), output_newhv] = fmincon(@(x_newhv) 0, x0, [], [], [], [], lb, ub, fun_newhv, options);
    % Extract solutions
    gamma_newhv(m) = x_newhv(1);
    I_sh_newhv(m) = x_newhv(2);
    I_r_newhv(m)=1/(Z_r+j*X_t)*(sqrt(3)*(V_s+V_se_newhv*exp(j*rho(m)))*exp(j*pi/6)/n_t-V_r*exp(j*theta_r));
    V_tap_newhv=0.5*(V_s+V_se_newhv*exp(j*rho(m)))*exp(j*-pi/3);
    Pr_newhv(m)=real(V_r*exp(j*theta_r)*conj( I_r_newhv(m)));
    Qr_newhv(m)=imag(V_r*exp(j*theta_r)*conj( I_r_newhv(m)));
    I_s_newhv(m)=0.5*I_sh_newhv(m)*exp(j*gamma_newhv(m))*exp(j*pi/3)+sqrt(3)/n_t* I_r_newhv(m)*exp(j*-pi/6);
    Ps_newhv(m)=real(V_s*conj(I_s_newhv(m)));
    Qs_newhv(m)=imag(V_s*conj(I_s_newhv(m)));
    Psh_newhv(m)=real(V_tap_newhv*conj(I_sh_newhv(m)*exp(j*gamma_newhv(m))));
    Qsh_newhv(m)=imag(V_tap_newhv*conj(I_sh_newhv(m)*exp(j*gamma_newhv(m))));
    Ssh_newhv(m) = sqrt(Psh_newhv(m)^2 + Qsh_newhv(m)^2);
    Pse_newhv(m)=real(V_se_newhv*exp(j*rho(m))*conj(I_sh_newhv(m)));
    Qse_newhv(m)=imag(V_se_newhv*exp(j*rho(m))*conj(I_sh_newhv(m)));
    Sse_newhv(m) = sqrt(Pse_newhv(m)^2 + Qse_newhv(m)^2);
    
end
polarplot( (Psh_newhv - Pse_newhv)/Srate,LineWidth=2)
hold on
polarplot( Qs_newhv/Srate,LineWidth=2)

% Nonlinear constraint function
function [c, ceq] = myfun_newhv(x_newhv, rho)
V_s = 10*1e3/sqrt(3); 
n_t = 25*sqrt(3);
V_r = 400/sqrt(3);
theta_r = pi/12;
Srate = 100*1e3/3;
Zbase = (V_r)^2/Srate;
X_r = 0.4*Zbase;   
X_t = 400^2/(300*1e3)*4/100;  
V_se_newhv = 0.2*V_s; 
Z_r = 1i * X_r;

% Extract variables
gamma_newhv = x_newhv(1);
I_sh_newhv = x_newhv(2);

% Compute intermediate quantities
I_r_newhv=1/(Z_r+j*X_t)*(sqrt(3)*(V_s+V_se_newhv*exp(j*rho))*exp(j*pi/6)/n_t-V_r*exp(j*theta_r));
I_s_newhv=0.5*I_sh_newhv*exp(j*gamma_newhv)*exp(j*pi/3)+sqrt(3)/n_t* I_r_newhv*exp(j*-pi/6);
Qs_newhv=imag(V_s*conj(I_s_newhv));
V_tap_newhv=0.5*(V_s+V_se_newhv*exp(j*rho))*exp(j*-pi/3);
Psh_newhv=real(V_tap_newhv*conj(I_sh_newhv*exp(j*gamma_newhv)));
Pse_newhv=real(V_se_newhv*exp(j*rho)*conj(I_s_newhv));
% Constraints
% Nonlinear equality constraints
ceq = [  Psh_newhv - Pse_newhv;Qs_newhv];
c = [];
end

Answer 1

[Torsten]

Use

% Options for fmincon
options = optimoptions('fmincon', 'Algorithm', 'interior-point', 'Display', 'off','ConstraintTolerance',1e-12);

instead of

% Options for fmincon
options = optimoptions('fmincon', 'Algorithm', 'interior-point', 'Display', 'off','ConstraintTolerance',1e-6);

But why do you use "fmincon" to solver a system of two nonlinear equations ? "fsolve" usually is the solver of choice.

rho = 0:2*pi/72:2*pi;
% Initial guess
x0 = [0 0];
% Options for fmincon
options = optimoptions('fmincon', 'Algorithm', 'interior-point', 'Display', 'off','ConstraintTolerance',1e-12);
% Constraints
lb = [];
ub = [];
% Solve for each rho
for m = 1:73
    % Define the nonlinear constraint function
    fun_newhv = @(x_newhv) myfun_newhv(x_newhv, rho(m));
    % Use fmincon to solve the problem
    [x_newhv, fval_newhv, exitflag_newhv(m), output_newhv] = fmincon(@(x_newhv) 0, x0, [], [], [], [], lb, ub, fun_newhv, options);
    % Extract solutions
    gamma_newhv(m) = x_newhv(1);
    I_sh_newhv(m) = x_newhv(2);
    [~,ceq(:,m)] = fun_newhv(x_newhv);
end
hold on
plot(1:m,ceq(1,:),'b')
plot(1:m,ceq(2,:),'r')
hold off
grid on

% Nonlinear constraint function
function [c, ceq] = myfun_newhv(x_newhv, rho)
V_s = 10*1e3/sqrt(3); 
n_t = 25*sqrt(3);
V_r = 400/sqrt(3);
theta_r = pi/12;
Srate = 100*1e3/3;
Zbase = (V_r)^2/Srate;
X_r = 0.4*Zbase;   
X_t = 400^2/(300*1e3)*4/100;  
V_se_newhv = 0.2*V_s; 
Z_r = 1i * X_r;
% Extract variables
gamma_newhv = x_newhv(1);
I_sh_newhv = x_newhv(2);
% Compute intermediate quantities
I_r_newhv=1/(Z_r+j*X_t)*(sqrt(3)*(V_s+V_se_newhv*exp(j*rho))*exp(j*pi/6)/n_t-V_r*exp(j*theta_r));
I_s_newhv=0.5*I_sh_newhv*exp(j*gamma_newhv)*exp(j*pi/3)+sqrt(3)/n_t* I_r_newhv*exp(j*-pi/6);
Qs_newhv=imag(V_s*conj(I_s_newhv));
V_tap_newhv=0.5*(V_s+V_se_newhv*exp(j*rho))*exp(j*-pi/3);
Psh_newhv=real(V_tap_newhv*conj(I_sh_newhv*exp(j*gamma_newhv)));
Pse_newhv=real(V_se_newhv*exp(j*rho)*conj(I_s_newhv));
% Constraints
% Nonlinear equality constraints
ceq = [  Psh_newhv - Pse_newhv;Qs_newhv];
c = [];
end    

[xin]

Thanks.The constraints are solved correctly.I mistyped expressions.I chose ''fmincon'' function because of the advice given in my previous question.Solver stopped prematurely. fsolve stopped because it exceeded the function evaluation limit, options.MaxFunctionEvaluations... - MATLAB Answers - MATLAB Central (mathworks.cn)

[Torsten]

rho = 0:2*pi/72:2*pi;
% Initial guess
x0 = [0 0];
% Options for fsolve
options = optimoptions('fsolve', 'Display', 'off');
% Solve for each rho
for m = 1:73
    % Define the equations
    fun_newhv = @(x_newhv) myfun_newhv(x_newhv, rho(m));
    % Use fsolve to solve the problem
    [x_newhv, fval_newhv, exitflag_newhv(m), output_newhv] = fsolve(fun_newhv, x0, options);
    % Extract solutions
    gamma_newhv(m) = x_newhv(1);
    I_sh_newhv(m) = x_newhv(2);
    res(:,m) = fval_newhv;
    x0 = x_newhv;
end
hold on
plot(1:m,res(1,:),'b')
plot(1:m,res(2,:),'r')
hold off
grid on

% Nonlinear equations
function res = myfun_newhv(x_newhv, rho)
V_s = 10*1e3/sqrt(3); 
n_t = 25*sqrt(3);
V_r = 400/sqrt(3);
theta_r = pi/12;
Srate = 100*1e3/3;
Zbase = (V_r)^2/Srate;
X_r = 0.4*Zbase;   
X_t = 400^2/(300*1e3)*4/100;  
V_se_newhv = 0.2*V_s; 
Z_r = 1i * X_r;
% Extract variables
gamma_newhv = x_newhv(1);
I_sh_newhv = x_newhv(2);
% Compute intermediate quantities
I_r_newhv=1/(Z_r+j*X_t)*(sqrt(3)*(V_s+V_se_newhv*exp(j*rho))*exp(j*pi/6)/n_t-V_r*exp(j*theta_r));
I_s_newhv=0.5*I_sh_newhv*exp(j*gamma_newhv)*exp(j*pi/3)+sqrt(3)/n_t* I_r_newhv*exp(j*-pi/6);
Qs_newhv=imag(V_s*conj(I_s_newhv));
V_tap_newhv=0.5*(V_s+V_se_newhv*exp(j*rho))*exp(j*-pi/3);
Psh_newhv=real(V_tap_newhv*conj(I_sh_newhv*exp(j*gamma_newhv)));
Pse_newhv=real(V_se_newhv*exp(j*rho)*conj(I_s_newhv));
res = [  Psh_newhv - Pse_newhv;Qs_newhv];
end