How can I find a unique solution within tolerance using solve?

Hello,

So I am using solve to find values for 3 variables while having 3 equations. The problem is that my equations dont seem to have unique solution. How can I setup tolerance in a way that the solution only exits within 20k to 2M? For example for values in 700k,1.5M and 2M the right hand side is approximately equal to the left hand side.

syms Rx Risop Rison
V1=252.79315;
V2=287.20685;
V3=496.51688;
V4=43.483124;
V5=38.272892;
V6=501.72711;
R=100000000;
R1=68220;
R2=68220;
Vx1=17.206846;
Vx2=226.51688;
Vx3=231.72711;
equ1=(V1/Risop)+(V1/R)==(Vx1/Rx)+(V2/Rison)+(V2/R);
equ2= (V3/Risop)+(V3/R)==(V4/Rison)+(V4/R1)+(V4/R)-(Vx2/Rx);
equ3= (V5/Risop)+(V5/R)+(V5/R2)==(V6/Rison)+(Vx3/Rx)+(V6/R);
[Risop,Rison,Rx]=solve([equ1,equ2,equ3],[Risop,Rison,Rx]);
vpa(Risop)
vpa(Rison)
vpa(Rx)
Risop=1500000;
Rison=2000000;
Rx=700000;
LHS=(V1/Risop)+(V1/R);
 RHS=(Vx1/Rx)+(V2/Rison)+(V2/R);
 vpa(LHS)
 vpa(RHS)
 LHS=(V3/Risop)+(V3/R);
 RHS=(V4/Rison)+(V4/R1)+(V4/R)-(Vx2/Rx);
 vpa(LHS)
 vpa(RHS)
 LHS=(V5/Risop)+(V5/R)+(V5/R2);
 RHS=(V6/Rison)+(Vx3/Rx)+(V6/R);
 vpa(LHS)
 vpa(RHS)

[Alex Sha]

Hi, how about the result below:

risop: -378351.252158428
rx: 124254.506405549
rison: -355908.211178662

[Taiji]

Yes, but this is not the solution. Risop=1500000; Rison=2000000; Rx=700000; should be the real solution that exist

Answer 1

[Torsten]

As long as rank(A) equals 3 in the below code, your system of linear equations has a unique solution.

syms Rx Risop Rison
V1=252.79315;
V2=287.20685;
V3=496.51688;
V4=43.483124;
V5=38.272892;
V6=501.72711;
R=100000000;
R1=68220;
R2=68220;
Vx1=17.206846;
Vx2=226.51688;
Vx3=231.72711;
equ1=(V1*Risop)+(V1/R)==(Vx1*Rx)+(V2*Rison)+(V2/R);
equ2= (V3*Risop)+(V3/R)==(V4*Rison)+(V4/R1)+(V4/R)-(Vx2*Rx);
equ3= (V5*Risop)+(V5/R)+(V5/R2)==(V6*Rison)+(Vx3*Rx)+(V6/R);
[A,b] = equationsToMatrix([equ1,equ2,equ3],[Risop,Rison,Rx]);
rank(A)

ans = 3

sol = solve([equ1,equ2,equ3],[Risop,Rison,Rx])

sol = struct with fields:
    Risop: -146640518737105837454502057739637344763690427/55481623593394354831199247574089475709215738691584
    Rison: -116915588320917442960783211679326820866617765/41611217695045766123399435680567106781911804018688
       Rx: 223257994367796429870241925302195019289661533/27740811796697177415599623787044737854607869345792

Risop = 1/sol.Risop;
Rison = 1/sol.Rison;
Rx = 1/sol.Rx;
format long
double([Risop;Rison;Rx])

ans = 3x1

1.0e+05 * -3.783512501947752 -3.559082094410595 1.242545059819753

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Taiji]

syms Rx Risop Rison positive
V1=252.79315;
V2=287.20685;
V3=496.51688;
V4=43.483124;
V5=38.272892;
V6=501.72711;
R=100000000;
R1=68220;
R2=68220;
Vx1=17.206846;
Vx2=226.51688;
Vx3=231.72711;
equ1=(V1*Risop)+(V1/R)==(Vx1*Rx)+(V2*Rison)+(V2/R);
equ2= (V3*Risop)+(V3/R)==(V4*Rison)+(V4/R1)+(V4/R)-(Vx2*Rx);
equ3= (V5*Risop)+(V5/R)+(V5/R2)==(V6*Rison)+(Vx3*Rx)+(V6/R);
[A,b] = equationsToMatrix([equ1,equ2,equ3],[Risop,Rison,Rx]);
rank(A)

ans = 3

sol1 = solve([equ1,equ2,equ3],[Risop,Rison,Rx])

sol1 = struct with fields:
    Risop: [0x1 sym]
    Rison: [0x1 sym]
       Rx: [0x1 sym]

Risop = 1/sol1.Risop;

Warning: Solution does not exist because the system is inconsistent.

Rison = 1/sol1.Rison;

Warning: Solution does not exist because the system is inconsistent.

Rx = 1/sol1.Rx;

Warning: Solution does not exist because the system is inconsistent.

format long
double([Risop;Rison;Rx])

ans =

  3x0 empty double matrix

here the rank is 3 but the solution is empty why? I used the Risop,Rison and Rx as positive. I do not want negative number solutions.

[Torsten]

If the rank is 3, there only exists one solution for your linear system of equations. If this unique solution has negative components (as is the case for your system) and you restrict the solution you want to be positive, you get an empty result.

If you want the "nearest" positive solution, use "lsqnonneg".

syms Rx Risop Rison 
V1=252.79315;
V2=287.20685;
V3=496.51688;
V4=43.483124;
V5=38.272892;
V6=501.72711;
R=100000000;
R1=68220;
R2=68220;
Vx1=17.206846;
Vx2=226.51688;
Vx3=231.72711;
equ1=(V1*Risop)+(V1/R)==(Vx1*Rx)+(V2*Rison)+(V2/R);
equ2= (V3*Risop)+(V3/R)==(V4*Rison)+(V4/R1)+(V4/R)-(Vx2*Rx);
equ3= (V5*Risop)+(V5/R)+(V5/R2)==(V6*Rison)+(Vx3*Rx)+(V6/R);
[A,b] = equationsToMatrix([equ1,equ2,equ3],[Risop,Rison,Rx]);
sol1 = lsqnonneg(double(A),double(b))

sol1 = 3x1

1.0e-05 * 0.1381 0.1214 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Risop = 1/sol1(1);
Rison = 1/sol1(2);
Rx = 1/sol1(3);
format long
double([Risop;Rison;Rx])

ans = 3x1

1.0e+05 * 7.241379477245834 8.235294079855887 Inf

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Answer 2

[Himanshu]

Hey Taiji,

To address the issue of non-unique solutions and ensure that the solutions fall within a specified range (20k to 2M), you can use MATLAB's numerical solver 'fsolve' from the Optimization Toolbox. The 'fsolve' function can be used to solve systems of nonlinear equations numerically and can handle bounds on the variables. By defining an initial guess and setting the lower and upper bounds for the variables, you can restrict the solutions to your desired range. Additionally, you can set tolerances for the function values and variables to ensure the solver's precision.

The code involves defining the equations as a function handle and specifying the initial guess for the variables. Then, we set the options for 'fsolve', including the tolerances and bounds. The 'fsolve' function is used to solve the system of equations, and the solutions are displayed along with a verification step to ensure the accuracy of the results. This approach allows you to handle the non-unique solutions and obtain results within the desired range.

You may go through the following documentaion link to know more about 'fsolve'.

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

% Define the variables and constants
V1 = 252.79315;
V2 = 287.20685;
V3 = 496.51688;
V4 = 43.483124;
V5 = 38.272892;
V6 = 501.72711;
R = 100000000;
R1 = 68220;
R2 = 68220;
Vx1 = 17.206846;
Vx2 = 226.51688;
Vx3 = 231.72711;

% Define the equations as functions
fun = @(x) [
    (V1/x(1)) + (V1/R) - (Vx1/x(3)) - (V2/x(2)) - (V2/R);
    (V3/x(1)) + (V3/R) - (V4/x(2)) - (V4/R1) - (V4/R) + (Vx2/x(3));
    (V5/x(1)) + (V5/R) + (V5/R2) - (V6/x(2)) - (Vx3/x(3)) - (V6/R)
];

% Initial guess for the variables
x0 = [1.5e6, 2e6, 700e3];

% Set options for fsolve, including bounds
options = optimoptions('fsolve', 'Display', 'iter', 'TolFun', 1e-6, 'TolX', 1e-6);
lb = [20e3, 20e3, 20e3]; % Lower bounds
ub = [2e6, 2e6, 2e6];   % Upper bounds

% Solve the system of equations
[x, fval, exitflag] = fsolve(fun, x0, options);

% Display the results
disp('Solutions:');
disp(['Risop = ', num2str(x(1))]);
disp(['Rison = ', num2str(x(2))]);
disp(['Rx = ', num2str(x(3))]);

% Verify the solutions
LHS1 = (V1/x(1)) + (V1/R);
RHS1 = (Vx1/x(3)) + (V2/x(2)) + (V2/R);
disp(['Equation 1: LHS = ', num2str(LHS1), ', RHS = ', num2str(RHS1)]);

LHS2 = (V3/x(1)) + (V3/R);
RHS2 = (V4/x(2)) + (V4/R1) + (V4/R) - (Vx2/x(3));
disp(['Equation 2: LHS = ', num2str(LHS2), ', RHS = ', num2str(RHS2)]);

LHS3 = (V5/x(1)) + (V5/R) + (V5/R2);
RHS3 = (V6/x(2)) + (Vx3/x(3)) + (V6/R);
disp(['Equation 3: LHS = ', num2str(LHS3), ', RHS = ', num2str(RHS3)]);

[John D'Errico]

@Himanshu

I'm sorry, but I think you don't understand how these tools work.

You defined the variables lb and ub. Good for you. But then you NEVER used them!

lb = [20e3, 20e3, 20e3]; % Lower bounds
ub = [2e6, 2e6, 2e6];   % Upper bounds
% Solve the system of equations
[x, fval, exitflag] = fsolve(fun, x0, options);

Just because you define variables called lb, and ub, does not mean that MATLAB will know what your intent is with those variables.

What is worse, is that fsolve does not offer bound constraints!

Instead, you could have suggested using lsqnonlin. it DOES allow bound constraints on the variables.