[python] Re: Stability Theory - bug :-(

At 08:45 PM 2/23/2005 +0100, you wrote:
good news 

I'm glad you think so!

 the results of your cool app now appear quite
realistic!


Even better...

The airbike data showed a stable range from 28 to 56 km/h
and P3 data showed no stable range at all.

I'll be gone for a couple of days (client trip), but when I get back I'll sort out the factors which contribute to the "B" and "E" instabilities. That should give an idea what/if they can be made self-stable. I really want to add in the human input factor as well as that pivot spring factor.

The output to stdout shows pairs of velocity and a negative
or positive number. The airbike output numbers are positive
from 17.5 to 34.5 km/h, then the output stops.
Why then is the stable range from 28 to 56 km/h?


The DOS output is still in MPH, sorry, I forgot to convert.

Three more features would be nice:
- input of data by file as a possible option


I was thinking of a standard .cfg type file last night

- different wheel sizes (don´t know if that matters at all)

no, only the momentum. It might be good to have wheel weight and diameter (assuming 95% is in the tire-rim) as the inputs, and calculate momentum internally instead.

- airbike has 20 cm of positive trail, GUI max are 10 cm 

! I didn't realize. I'll change it. Any other limit problems?

Have you thought of making the sourcecode available
for everyone?

Yes, I often do Python (the language) stuff LGPL, only reserving code authorship. The algorithm below follows Hand's. If anyone can or wants to program in Python, they could easily extend it. I use the Boa Constructor IDE.
I've been working on a larger cycling software package of my own, so I plan on including a version of this as a module, and so part of the impetus for the graphics and GUI. If anyone wants to collaborate, I'd be quite happy, and anything anyone adds would of course then belong to them and/or the group. It's a lot more fun and rewarding than writing boring apps for corporate clients.

The main algorithm is RouthHurwitz (the simpler DOS output is included, current code is at home).
In this version, lines 62-76 iterate over desired input ranges, you could make them singular. All are imperial units, and it writes to a file; the print lines are just commented out 165-172.
If you trace back the b0 and e0 formulas (line 130+), you'll see what factors are driving them negative, which is the criteria defining Routh-Hurwitz instability.

Cheers,
Ray

#-----------------------------------------------------------------------------
# Name:        FrameGeo.py
# Purpose:
#
# Author:      Ray Schumacher
#
# Created:     2005/23/02
# RCS-ID:      $Id: FrameGeo.py $
# Copyright:   (c) 2005
# Licence:     Not-for-profit public use
#-----------------------------------------------------------------------------
#Boa:Frame:wxFrameGeometry

import math
import time

def __test():

    outFile=open('results.csv',"w")
    outFile.write(  'rearMassTrail'+','+
                    'htRearMass'+','+
                    'fM2Piv'+','+
                    #'frtMassHt'+','+
                    'wheelbase'+','+
                    'headAngle'+','+
                    'forkTrail'+','+
                    'stableRange[0]'+','+
                    'stableRange[-1]'+','+
                    'stableRange[-1]-stableRange[0]'+'\n')
    rearMass = 170.0    ## lbm
    rearMassTrail = 18. ## in
    htRearMass = 14.0   ## in
    Ryp = 100000.0      ## lbm in**2
    Rzp = 20000.0       ## lbm in**2
    alphar = 60.0 * math.pi / 180.0    ## degress to radians
    wheelbase = 52.0    ## in
    forkTrail = -2.0    ## in
    headAngle = 35.0 * math.pi / 180.0 ## degress to radians
    frontMass = 30.0    ## lbm
    fM2Piv = 12.0       ## in
    frtMassHt = 14. ## in
    Fyp = 50.0      ## lbm in**2
    Fzp = 2000.0    ## lbm in**2
    alphaf = 45.0 * math.pi / 180.0    ## degress to radians
    CrR = 50.0      ## lbm in
    CfR = 50.0      ## lbm in
    g = 386.4

    ## inertias
    Ryy = Rzp*(math.cos(alphar))**2+Ryp*(math.sin(alphar))**2
    Rzz = Rzp*(math.sin(alphar))**2+Ryp*(math.cos(alphar))**2
    Ryz = Rzp*math.cos(alphar)*math.sin(alphar)- \
            Ryp*math.sin(alphar)*math.cos(alphar)

    Fyy = Fzp*(math.cos(alphaf))**2+Fyp*(math.sin(alphaf))**2
    Fzz = Fzp*(math.sin(alphaf))**2+Fyp*(math.cos(alphaf))**2
    Fyz = Fzp*math.cos(alphaf)*math.sin(alphaf)- \
            Fyp*math.sin(alphaf)*math.cos(alphaf)

    startTime = time.time()

    for thisRMassTr in range(22, 27):
        rearMassTrail = float(thisRMassTr)
        print '\n', rearMassTrail
        for thishtRearMass in range(15, 20):
            htRearMass = float(thishtRearMass)
            for thisfM2Piv in range(5, 15):
                fM2Piv = float(thisfM2Piv)
                for thisfrtMassHt in range(14, 18):
                    frtMassHt = float(thisfrtMassHt)
                    lastTime = time.time()
                    for thiswheelbase in range(50, 60, 2):
                        wheelbase = float(thiswheelbase)
                        for thisAngle in range(27, 65):
                            headAngle = float(thisAngle)
                            for thisTrail in range(-12, -6):
                                stableRange = []
                                forkTrail = thisTrail/2.

                                ## front lengths
                                hf = 
(forkTrail+fM2Piv)*math.sin(headAngle)+frtMassHt*math.cos(headAngle)
                                lf = 
(forkTrail+fM2Piv)*math.cos(headAngle)-frtMassHt*math.sin(headAngle)

                                ## total bike inertias and special front inertia
                                Tyy = 
rearMass*htRearMass**2+Ryy+frontMass*hf**2+Fyy
                                Tzz = 
rearMass*rearMassTrail**2+Rzz+frontMass*(wheelbase+lf)**2+Fzz
                                Tyz = 
-rearMass*htRearMass*rearMassTrail+Ryz-frontMass*hf*(wheelbase+lf)+Fyz

                                Fll = 
frontMass*fM2Piv**2+Fyy*(math.sin(headAngle))**2 - \
                                        
Fyz*math.sin(2*headAngle)+Fzz*(math.cos(headAngle))**2
                                Fly = 
-frontMass*hf*fM2Piv-Fyy*math.sin(headAngle)+Fyz*math.cos(headAngle)
                                Flz = 
frontMass*(wheelbase+lf)*fM2Piv-Fyz*math.sin(headAngle)+Fzz*math.cos(headAngle)

                                ## design length ratios
                                Ktr = forkTrail/wheelbase
                                Khd = math.cos(headAngle)/wheelbase
                                mt = rearMass+frontMass
                                ht = (rearMass*htRearMass+frontMass*hf)/mt
                                lt = 
(rearMass*rearMassTrail+frontMass*(wheelbase+lf))/mt
                                nu = frontMass*fM2Piv+mt*lt*Ktr

                                """print Ryy, Rzz,Ryz
                                print Fyy,Fzz , Fyz
                                print Tyy,Tzz,Tyz
                                print Fll,Fly,Flz
                                print Ktr,Khd,mt
                                print ht,lt,nu"""

                                ## coeffs for the lean equation
                                a1 = Tyy
                                a2 = -g*mt*ht
                                a3 = Fly+Ktr*Tyz
                                a4 = 
-(CfR*math.cos(headAngle)+Ktr*(CfR+CrR))+(Tyz*Khd-Ktr*mt*ht)
                                a5 = g*nu
                                a6 = -(CfR+CrR)*Khd-Khd*mt*ht
                                #print a1,a2,a3 ,a4 ,a5 , a6

                                ## coeffs for the steer equation
                                b1 = Fll+2*Ktr*Flz+(Ktr**2)*Tzz
                                b2 = Khd*(Flz+Ktr*Tzz)+Ktr*nu
                                b3 = -g*nu*math.sin(headAngle)
                                b4 = CfR*math.sin(headAngle)*Khd+nu*Khd
                                b5 = Fly+Ktr*Tyz
                                b6 = CfR*math.cos(headAngle)+Ktr*(CfR+CrR)
                                b7 = g*nu
                                #print b1,b2,b3,b4,b5,b6,b7

                                ## coeffs for the 4th order polynomial
                                a0 = a1*b1-a3*b5
                                b0 = a1*b2-a3*b6-a4*b5
                                c0 = a1*b3-a3*b7-a5*b5+b1*a2
                                d0 = a1*b4-b5*a6-b6*a4
                                e0 = -a4*b7-a5*b6+b2*a2
                                f0 = -b6*a6
                                g0 = a2*b3-a5*b7
                                h0 = a2*b4-a6*b7
                                i0 = b0*d0*f0-a0*f0**2
                                j0 = b0*c0*f0+b0*d0*e0-2*a0*e0*f0-h0*b0**2
                                k0 = b0*c0*e0-a0*e0**2-g0*b0**2
                                l0 = b0*d0-a0*f0
                                m0 = b0*c0-a0*e0

                                """print 'A =',a0,'>0'
                                print 'B =',b0,'V>0'
                                print 'C =',c0,'+',d0,'V**2>0'
                                print 'D =',e0,'V+',f0, 'V**3>0'
                                print 'E =',g0,'+',h0,'V**2>0'
                                print '6th =',l0,'V**3+',m0, 'V>0'
                                print '7th =',i0,'V**4+',j0,'V**2+', k0, '>0'"""


                                ## print stability
                                for j in range(0, 100):
                                    PRT = [' ']*7
                                    speed = j
                                    Vmph = speed/2.
                                    Vips = Vmph*5280*12/3600.
                                    A = a0
                                    B = b0*Vips
                                    C = c0+d0*Vips**2
                                    bigD = e0*Vips+f0*Vips**3
                                    E = g0+h0*Vips**2
                                    rOuth6 = l0*Vips**3 +m0*Vips
                                    rOuth7 = i0*Vips**4+j0*Vips**2+k0
                                    """if A<=0.: PRT[0] = 'A'
                                    if B<=0.: PRT[1] = 'B'
                                    if C<=0.: PRT[2] = 'C'
                                    if bigD<=0.: PRT[3] = 'D'
                                    if E<=0.: PRT[4] = 'E'
                                    if rOuth6<=0.: PRT[5] = '6'
                                    if rOuth7<=0.: PRT[6] = '7'
                                    print 'Vel=', round(Vmph, 2), PRT[0], 
PRT[1], PRT[2], PRT[3], PRT[4], PRT[5], PRT[6], rOuth7"""

                                    if min([A,B,C,bigD,E,rOuth6,rOuth7])>0:
                                        stableRange.append(Vmph)


                                if len(stableRange) and 
stableRange[-1]-stableRange[0]>10 and stableRange[0]<10:
                                    outFile.write(str(rearMassTrail)+','+
                                        str(htRearMass)+','+
                                        str(fM2Piv)+','+
                                        #str(frtMassHt)+','+
                                        str(wheelbase)+','+
                                        str(headAngle)+','+
                                        str(forkTrail)+','+
                                        str(stableRange[0])+','+
                                        str(stableRange[-1])+','+
                                        
str(stableRange[-1]-stableRange[0])+'\n')
                                    #print '.',
                                    #print 
rearMassTrail,htRearMass,fM2Piv,frtMassHt,wheelbase,headAngle,forkTrail, '\t', 
stableRange[0], '-', stableRange[-1], '\t', stableRange[-1]-stableRange[0]
                                    #print headAngle, forkTrail, '\t', 
stableRange[0], '-', stableRange[-1], '\t', stableRange[-1]-stableRange[0]

                    print round((97*(time.time()-lastTime) - 
(time.time()-startTime))/60, 1), 'mins remain'
    outFile.close()

if __name__ == '__main__':
    __test()

Other related posts: