Solver stopped prematurely. fsolve stopped because it exceeded the function evaluation limit, options.Ma​xFunctionE​valuations = 2.000000e+02.

I want to solve a system of equations containing two non-linear equations, one of which is Psh==Pse,the other is Qs==0, but solving it with the fsolve function doesn't achieve the given constraints, how can I solve this problem?

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/(150*1e3)*4/100;

rho=0:2*pi/24:2*pi;

V_se=0.2*V_r;

Z_r=j*X_r;

for m=1:25

fun=@(x)myfun(x,rho(m));

x0=[0 0];

[x,fval,exitflag(m),~]=fsolve(fun,x0);

gamma(m)=x(1);

I_sh(m)=x(2);

I_r(m)=1/(Z_r+j*X_t)*(sqrt(3)*V_s*exp(j*pi/6)/n_t+V_se*exp(j*rho(m))-V_r*exp(j*theta_r)-j*X_t*I_sh(m)*exp(j*gamma(m)));

V_c1(m)=sqrt(3)*V_s*exp(j*pi/6)/n_t-j*X_t*(I_sh(m)*exp(j*gamma(m))+I_r(m));

Psh(m)=real(V_c1(m)*conj(I_sh(m)*exp(j*gamma(m))));

Pse(m)=real(V_se*exp(j*rho(m))*conj(I_r(m)));

Pr(m)=real(V_r*exp(j*theta_r)*conj(I_r(m)));

Qr(m)=imag(V_r*exp(j*theta_r)*conj(I_r(m)));

Sr(m)=sqrt(Pr(m)^2+Qr(m)^2);

Ps(m)=real(sqrt(3)*V_s*exp(j*pi/6)/n_t*conj(I_r(m)+I_sh(m)*exp(j*gamma(m))));

Qs(m)=imag(sqrt(3)*V_s*exp(j*pi/6)/n_t*conj(I_r(m)+I_sh(m)*exp(j*gamma(m))));

Ss(m)=sqrt(Ps(m)^2+Qs(m)^2);

end

plot( Ps/Srate,Qs/Srate,LineWidth=2,color='[0.9290 0.6940 0.1250]')

function y=myfun(x,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/(150*1e3)*4/100;

V_se=0.2*V_r;

Z_r=j*X_r;

%%%%%%%%%%%%%%%%Two constraints,Psh==Pse,Qs==0

gamma=x(1);

I_sh=x(2);

I_r=1/(Z_r+j*X_t)*(sqrt(3)*V_s*exp(j*pi/6)/n_t+V_se*exp(j*rho)-V_r*exp(j*theta_r)-j*X_t*I_sh*exp(j*gamma));

V_c1=sqrt(3)*V_s*exp(j*pi/6)/n_t-j*X_t*(I_sh*exp(j*gamma)+I_r);

Psh=real(V_c1*conj(I_sh*exp(j*gamma)));

Pse=real(V_se*exp(j*rho)*conj(I_r));

Qs=imag(sqrt(3)*V_s*exp(j*pi/6)/n_t*conj(I_r+I_sh*exp(j*gamma)));

y=[Psh-Pse; Qs];

end

Answer 1

[Anagha Mittal]

Hi,

"fsolve" function is not giving you the desired solution as it is not able to handle the constraints. I would rather suggest to use "fmincon" function. I have made a few changes in your code accordingly, please take a look below:

% 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/(150*1e3)*4/100;  
rho = 0:2*pi/24:2*pi; 
V_se = 0.2*V_r; 
Z_r = 1i * X_r;

% Initial guess
x0 = [0 0];

% Options for fmincon
options = optimoptions('fmincon', 'Algorithm', 'interior-point', 'Display', 'iter');

% Constraints
lb = [];
ub = [];

% Solve for each rho
for m = 1:25
    % Define the nonlinear constraint function
    fun = @(x) myfun(x, rho(m));
    
    % Use fmincon to solve the problem
    [x, fval, exitflag(m), output] = fmincon(@(x) 0, x0, [], [], [], [], lb, ub, fun, options);
    
    % Extract solutions
    gamma(m) = x(1);
    I_sh(m) = x(2);
    
    % Compute derived quantities
    I_r(m) = 1 / (Z_r + 1i * X_t) * (sqrt(3) * V_s * exp(1i * pi/6) / n_t + V_se * exp(1i * rho(m)) - V_r * exp(1i * theta_r) - 1i * X_t * I_sh(m) * exp(1i * gamma(m)));
    V_c1(m) = sqrt(3) * V_s * exp(1i * pi/6) / n_t - 1i * X_t * (I_sh(m) * exp(1i * gamma(m)) + I_r(m));
    Psh(m) = real(V_c1(m) * conj(I_sh(m) * exp(1i * gamma(m))));
    Pse(m) = real(V_se * exp(1i * rho(m)) * conj(I_r(m)));
    Pr(m) = real(V_r * exp(1i * theta_r) * conj(I_r(m)));
    Qr(m) = imag(V_r * exp(1i * theta_r) * conj(I_r(m)));
    Sr(m) = sqrt(Pr(m)^2 + Qr(m)^2);
    Ps(m) = real(sqrt(3) * V_s * exp(1i * pi/6) / n_t * conj(I_r(m) + I_sh(m) * exp(1i * gamma(m))));
    Qs(m) = imag(sqrt(3) * V_s * exp(1i * pi/6) / n_t * conj(I_r(m) + I_sh(m) * exp(1i * gamma(m))));
    Ss(m) = sqrt(Ps(m)^2 + Qs(m)^2);
end

% Plot results
plot(Ps/Srate, Qs/Srate, 'LineWidth', 2, 'Color', [0.9290 0.6940 0.1250]);

% Nonlinear constraint function
function [c, ceq] = myfun(x, 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/(150*1e3)*4/100;  
    V_se = 0.2*V_r; 
    Z_r = 1i * X_r;
    
    % Extract variables
    gamma = x(1);
    I_sh = x(2);
    
    % Compute intermediate quantities
    I_r = 1 / (Z_r + 1i * X_t) * (sqrt(3) * V_s * exp(1i * pi/6) / n_t + V_se * exp(1i * rho) - V_r * exp(1i * theta_r) - 1i * X_t * I_sh * exp(1i * gamma));
    V_c1 = sqrt(3) * V_s * exp(1i * pi/6) / n_t - 1i * X_t * (I_sh * exp(1i * gamma) + I_r);
    
    % Constraints
    Psh = real(V_c1 * conj(I_sh * exp(1i * gamma)));
    Pse = real(V_se * exp(1i * rho) * conj(I_r));
    Qs = imag(sqrt(3) * V_s * exp(1i * pi/6) / n_t * conj(I_r + I_sh * exp(1i * gamma)));
    
    % Nonlinear equality constraints
    ceq = [Psh - Pse; Qs];
    c = [];
end

For more information on "fmincon", refer to the following documentation:

https://www.mathworks.com/help/optim/ug/fmincon.html

Also, you may refer to the following documentation to read about constraints while solving non-linear equations:

https://www.mathworks.com/help/optim/ug/nonlinear-systems-with-constraints.html

Hope this helps!

[xin]

Thanks，you really help me a lot and the code is written cleanly and easy to understand.