Python Forum
Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
Solve system of equations
#11
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.
Reply
#12
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.
Reply
#13
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.
Reply
#14
Quote: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.
Reply
#15
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?
Reply
#16
(Oct-26-2020, 10:05 AM)Askic Wrote: 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
Reply
#17
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.
Reply
#18
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)
Reply
#19
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???
Reply
#20
(Oct-26-2020, 05:28 PM)Askic Wrote: 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.
Gribouillis likes this post
Reply


Possibly Related Threads…
Thread Author Replies Views Last Post
Heart how to solve complex equations in python HoangF 3 2,731 Dec-26-2021, 07:04 PM
Last Post: HoangF
Question Help with code to solve polynomials equations hiviera 1 1,761 Jul-31-2021, 01:56 AM
Last Post: hiviera
  Difference between os.system("clear") and os.system("cls") chmsrohit 7 16,492 Jan-11-2021, 06:30 PM
Last Post: ykumar34
  How to solve equations, with groups of variables and or constraints? ThemePark 0 1,644 Oct-05-2020, 07:22 PM
Last Post: ThemePark
  Solve a system of linear equations with binary variables lopeslimagabriel 3 2,443 Sep-24-2020, 07:09 AM
Last Post: scidam
  How to solve difficult non-linear equations? shreeniket 3 2,348 Apr-23-2020, 01:36 AM
Last Post: shreeniket
Question Difference between Python's os.system and Perl's system command Agile741 13 6,650 Dec-02-2019, 04:41 PM
Last Post: Agile741
  System of 3 non-linear equations in 3 unknowns (how-to-solve?) samsonite 2 3,537 Mar-23-2019, 10:14 AM
Last Post: samsonite
  Beginner: System of Equations Mahdi1994 2 2,465 Mar-19-2018, 12:37 AM
Last Post: Tiskolin

Forum Jump:

User Panel Messages

Announcements
Announcement #1 8/1/2020
Announcement #2 8/2/2020
Announcement #3 8/6/2020