Solve system of equations

[Sancho_Pansa]

I've tried with fsolve.
Here is my code:

import math
import numpy as np
from scipy.optimize import fsolve

def f(LCR):
    #variables
    L  = LCR[0]
    Cs = LCR[1]
    Rs = LCR[2]
    
    # global constants
    Pi = math.pi
    Rin = 10.0
    Cl = 3.0e-12
    
    # measurement-related constants
    f0 = 130e6
    f1 = 140e6
    f2 = 150e6
    omg = [None]*3
    out = [None]*3
    omg[0] = 2*Pi*f0
    omg[1] = 2*Pi*f1
    omg[2] = 2*Pi*f2
    out[0] = 6.658482971121064
    out[1] = 25.89226207742419
    out[2] = 7.013950018010094

    # A and B arrays
    Ar = [None]*3
    Aj = [None]*3
    Br = [None]*3
    Bj = [None]*3
    equ = [None]*3
    
    # equation building
    for i in range(3):
        Ar[i] = 1 - omg[i]**2*Cl*L
        Aj[i] = omg[i]*Cs*Rs
        Br[i] = 1 - omg[i]**2*(L*(Cs + Cl) + Cl*Cs*Rin*Rs)
        Bj[i] = omg[i]*(Cs*Rs + Cl*Rin + Cl*Rs - omg[i]**2*Cl*Cs*Rin*L)

        equ[i] = out[i] - math.sqrt(Ar[i]*Ar[i] + Aj[i]*Aj[i])/math.sqrt(Br[i]*Br[i] + Bj[i]*Bj[i])

    return equ

LCR0 = np.array([0.3e-6, 1e-12, 1])
LCR = fsolve(f, LCR0)
LCR

Here is output:

Output:
Warning (from warnings module):
  File "/home/pavel47/.local/lib/python3.8/site-packages/scipy/optimize/minpack.py", line 175
    warnings.warn(msg, RuntimeWarning)
RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last five Jacobian evaluations.
>>> LCR
array([ 3.10756682e-07,  1.02921301e-12, -4.22520679e+00])
>>> 

While the first two values are correct, the third is completely wrong.
Any suggestions ?
Thanks.

[Sancho_Pansa]

Here is updated version of code:

import math
import numpy as np
from scipy.optimize import fsolve

def f(LCR, *args):
    omg = np.array([args[0], args[1], args[2]])
    out = np.array([args[3], args[4], args[5]])
   
    # variables
    L  = LCR[0]
    Cs = LCR[1]
    Rs = LCR[2]
    
    # global constants
    Rin = 10.0
    Cl = 3.0e-12

    # declaration A, B and equ arrays
    Ar, Aj, Br, Bj, equ = [np.zeros(3) for _ in range(5)]
    
    # equation building
    for i in range(3):
        Ar[i] = 1 - omg[i]**2*Cl*L
        Aj[i] = omg[i]*Cs*Rs
        Br[i] = 1 - omg[i]**2*(L*(Cs + Cl) + Cl*Cs*Rin*Rs)
        Bj[i] = omg[i]*(Cs*Rs + Cl*Rin + Cl*Rs - omg[i]**2*Cl*Cs*Rin*L)

        #equ[i] = out[i] - math.sqrt(Ar[i]*Ar[i] + Aj[i]*Aj[i])/math.sqrt(Br[i]*Br[i] + Bj[i]*Bj[i])
        equ[i] = out[i]*math.sqrt(Br[i]*Br[i] + Bj[i]*Bj[i]) - math.sqrt(Ar[i]*Ar[i] + Aj[i]*Aj[i])
    return equ

# Define frequencies in MHz
f1 = 130
f2 = 140
f3 = 150
omega = np.array([f1, f2, f3])*1.0e6*2*math.pi

# Define outputs:
out1 = 6.658482971121064
out2 = 25.89226207742419
out3 = 7.013950018010094
out = np.array([out1, out2, out3])

# Define guesses for L [uH], Cs [pF], Rs [Ohm]
L0  = 0.3e-6
Cs0 = 1e-12
Rs0 = 5

LCR0 = np.array([L0, Cs0, Rs0])
LCR = fsolve(f, LCR0, args=tuple(np.concatenate((omega, out))), maxfev=10000)
print(LCR)

With initial conditions closer to the solution, the result is much worse.

[Gribouillis]

Unfortunately, I don't have much time to help you now but I think you reduced too fast the complex equations to real numbers. Go back to the complex equations and try to interprete them in geometrical terms. It seems to me that with some work you could express Cs in terms of L and Rs, then design an optimisation form of the problem. When I have the time to do so, I may answer in more details but I'm afraid it won't be before a week because I'm going to be on a journey for a few days.

[Sancho_Pansa]

It seems to me that with some work you could express Cs in terms of L and Rs

Frankly speaking I didn't understand your idea. Well ... after hard work I could do that. And after that ... what?
I think a solution can come from some parameters passing in fsolve.

[Askic]

Sancho_Pansa,

can you please elaborate more about the problem itself. You said you have 3 equations with 3 unknowns. The unknown variables that need to be found are: Cs, L, Rs
So the first equation is the equation for Br, second is for Bj and the third equation is for TF.
When you say that frequency omega is different, you mean in the first equation is omega1, in second is omega2 and in the third one is omega3 or maybe something else?

[Sancho_Pansa]

Sancho_Pansa,

can you please elaborate more about the problem itself. You said you have 3 equations with 3 unknowns. The unknown variables that need to be found are: Cs, L, Rs
So the first equation is the equation for Br, second is for Bj and the third equation is for TF.
When you say that frequency omega is different, you mean in the first equation is omega1, in second is omega2 and in the third one is omega3 or maybe something else?

Well, at the beginning I badly formulated the problem.
Here is correct setup:

• there are 3 variables: Cs, L, Rs
  
• there are 3 equations: How these 3 equations differ - at each frequency there is a specific output
  

[Askic]

I see,
this looks like as if the constrained optimization would probably do better (since L,C and R cannot be negative)

import math
from scipy.optimize import least_squares


def f(LCR):
    # variables
    L = LCR[0]
    Cs = LCR[1]
    Rs = LCR[2]

    # global constants
    Pi = math.pi
    Rin = 10.0
    Cl = 3.0e-12

    # measurement-related constants
    f0 = 130e6
    f1 = 140e6
    f2 = 150e6
    omg = [0, 0, 0]
    out = [0, 0, 0]
    omg[0] = 2*Pi*f0
    omg[1] = 2*Pi*f1
    omg[2] = 2*Pi*f2
    out[0] = 6.658482971121064
    out[1] = 25.89226207742419
    out[2] = 7.013950018010094

    # A and B arrays
    Ar = [0, 0, 0]
    Aj = [0, 0, 0]
    Br = [0, 0, 0]
    Bj = [0, 0, 0]
    equ = [0, 0, 0]

    # equation building
    for i in range(3):
        Ar[i] = 1 - omg[i]**2*Cl*L
        Aj[i] = omg[i]*Cs*Rs
        Br[i] = 1 - omg[i]**2*(L*(Cs + Cl) + Cl*Cs*Rin*Rs)
        Bj[i] = omg[i]*(Cs*Rs + Cl*Rin + Cl*Rs - omg[i]**2*Cl*Cs*Rin*L)

        equ[i] = out[i] - math.sqrt(Ar[i]*Ar[i] + Aj[i]
                                    * Aj[i])/math.sqrt(Br[i]*Br[i] + Bj[i]*Bj[i])
    return equ


LCR0 = [0.3e-6, 1e-12, 1]
res = least_squares(f, (LCR0), bounds=((0, 0, 0), (1.0e-6, 1.0e-10, 2)))
print(res)

Let me know if this is closer to the expected solution.

[Sancho_Pansa]

Thanks.
1st value (i.e. L) is Ok.
2nd value (i.e. C) is lower than expected of 30%
The last one (i.e. R) is 40% of expected.
Here is code, that generates expected values:

import math

Pi = math.pi
for f_raw in (130, 140, 150):
    Cs_raw = 1.0
    Cl_raw = 3.0
    L_raw = 0.32

    Cs = Cs_raw*1e-12
    Cl = Cl_raw*1e-12
    L = L_raw*1e-6
    Rs = 5.0
    Rin = 10.0
    f = f_raw*1e6
    omega = 2*Pi*f

    Ar = 1 - omega**2*Cs*Cl
    Aj = omega*Cs*Rs
    A = math.sqrt(Ar**2 + Aj**2)

    Br = 1 - omega**2*Cs*L - omega**2*Cl*L - omega**2*Cl*Cs*Rin*Rs
    Bj = omega*Cs*Rs + omega*Cl*Rin - omega**3*Cl*Cs*Rin*L + omega*Cl*Rs
    B = math.sqrt(Br**2 + Bj**2)

    print(f_raw, A/B)

[Askic]

Hello Sancho,

please double check your equations. I'm looking your last code. Is it Ar = 1 - omega**2*Cs*Cl or
Ar = 1 - omega**2*Cl*L???

[Sancho_Pansa]

Hello Sancho,

please double check your equations. I'm looking your last code. Is it Ar = 1 - omega**2*Cs*Cl or
Ar = 1 - omega**2*Cl*L???

Hello Askic.
Here is the solution, which works fine.
Modifications that was made:

• I replaced the capacitance Cl by the resistor Rl which drastically simplifies the transfer function, which in turn simplifies the work of the fsolve function.
  
• The equations are now generated by transformations performed on the circuit, created by the Circuit function, which minimizes typos.
  

from lcapy import Circuit, j, omega
import numpy as np
from math import pi
from scipy.optimize import fsolve

def f(LCR, *args):
    omg = np.array([args[0], args[1], args[2]])
    out = np.array([args[3], args[4], args[5]])

    cct = Circuit('''
    ... Rin 1 2
    ... Cs  2 4
    ... La  2 3
    ... Rs  3 4
    ... Rl  4 0''')

    cct1 = cct.subs({'Rin': 50, 'Cs':LCR[1], 'La': LCR[0], 'Rs': LCR[2], 'Rl': 100})

    H = cct1.transfer(1,0,4,0)
    A = H(j*omega).simplify()
    Am = A.magnitude

    equ = np.zeros(3)
    
    # equation building
    for i in range(3):
        equ[i] = out[i] - Am.evaluate(omg[i])
    return equ

# Define frequencies in MHz
f1 = 130
f2 = 140
f3 = 150
omg = np.array([f1, f2, f3])*1.0e6*2*pi

# Define outputs:
out = [0.27175539, 0.24606606, 0.22193977]

# Define guesses for L [uH], Cs [pF], Rs [Ohm]
L0  = 0.3e-6
Cs0 = 1e-12
Rs0 = 1

LCR0 = np.array([L0, Cs0, Rs0])
LCR = fsolve(f, LCR0, args=tuple(np.concatenate((omg, out))), maxfev=10000)
print(LCR)

Here is output:

Output:
[3.20000543e-07 9.99996130e-13 4.99951198e+00]

Here are real values:
0.32e-6, 1.0e-12, 5.0

Sincerely,

Sancho.