Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

Register Log in
  • Congratulations cowski on being selected by the Eng-Tips community for having the most helpful posts in the forums last week. Way to Go!

Python Solve Functions - Help For A Noobie 2

Status
Not open for further replies.
sticksandtriangles

sticksandtriangles

Structural
Joined
Apr 7, 2015
Messages
500
Location
US
sparked by the python coding discussion by Celt83, I am diving in and trying to solve instantaneous center of rotation problems with python.

The example problem I am trying to solve is outlined here:

The way I see this problem, there are (3) unknowns:
[ol 1]
[li]X location of the instantaneous center of rotation, [/li]
[/ol]
[ol 2]
[li]Y location of the instantaneous [/li]
[/ol]
[ol 3]
[li]and the Applied Loading that needs to be maxed out.[/li]
[/ol]

(sorry don't know how to fix the numbering above, fiddled for too damn long)

I also have (3) equations that need to equal 0
[ul]
[li]Sum of Forces in the X = 0[/li]
[/ul]
[ul]
[li] Sum of Forces in the Y = 0[/li]
[/ul]
[ul]
[li]Sum of Moments = 0 [/li]
[/ul]



Even in my "native coding language" matlab, this would be somewhat challenging to code, so I am sure it will be fun using this as a learning experience in python.

Anyways, python experts, do you have any suggestions on how to solve this? Should I explore symbolic notation within in python? Should I use a solve block that iterates through the (3) unknowns?

Any help or suggestions/starting points would be appreciated.

My base code that does the guess and check version of this is shown below:
image_one_cr58iv.png

The final print shows the comparison of applied loading verse resist force.
image_2_fuboii.png







S&T
 
Sort by date Sort by votes
The Scipy library has several solver functions that will do what you want.

Examples can be found in my xlScipy3 spreadsheet, which can be downloaded from:

The spreadsheet calls the Python functions from Excel using xlwings, but if you prefer to work directly with Python the Python code is included in the download.

For solvers with more than one unknown see the Solvers2 sheet, which use the Scipy.Optimize minimize function. You need to write a function that returns the difference between calculated and target values for given inputs, with output as a Numpy array, and the minimize function will adjust the inputs to minimize the differences.



Doug Jenkins
Interactive Design Services
 
Upvote 0 Downvote
IDS, thank you.

I am digging into these spreadsheets now and will check in later (got scipy downloaded).
Feel a little overwhelmed currently [spin2], lots of learning to ensue...

S&T
 
Upvote 0 Downvote
Let us know how you go.

To use Scipy you will also need Numpy, but they are both pretty well indispensable for engineering work with Python.

As an example of what you can do with Numpy, using the xloc list from your code:

Python: 
import numpy as np
…
# Convert xloc to a Numpy array
xloc = np.array(xloc)
ICX = -1.69
# Subtract ICX from all xloc values in one operation
xIC = xloc - ICX

Doug Jenkins
Interactive Design Services
 
Upvote 0 Downvote
hmm learned something new hadn't ever looked into Scipy for this, here's my code below which seems to check out:
It is pretty sensitive to where the initial IC coordinates are assumed to be

As IDS notes you can condense a lot of this down to single line steps using numpy arrays, I kept the bulk of it as native python lists and broken up into the calculation steps.

Python: 
from __future__ import division
import math as m
import scipy.optimize as sci

def ic_check(IC,xloc,yloc, ecc, theta):
    num_bolts = len(xloc)
    
    deltamax = 0.34
    Rult = 32.471
    
    ICx = IC[0]
    ICy = IC[1]
    xIC = []
    yIC = []
    di = []
    deltai = []
    ri = []
    
    fx = []
    fy = []
    moment = []
    
    for x in xloc:
        xICtemp = x - ICx
        xIC.append(xICtemp)
        
    for y in yloc:
        yICtemp = y - ICy
        yIC.append(yICtemp)
    
    i=0
    for i in range(num_bolts):
        ditemp = m.sqrt((xIC[i]*xIC[i])+(yIC[i]*yIC[i]))
        di.append(ditemp)
    
    dmax = max(di)
    
    i=0
    for i in range(num_bolts):
        deltaitemp = (di[i]/dmax)*deltamax
        deltai.append(deltaitemp)
        
    i=0
    for i in range(num_bolts):
        ritemp = Rult*m.pow(1-m.pow(m.e,-10.0*deltai[i]),0.55)
        ri.append(ritemp)
        
    i=0
    for i in range(num_bolts):    
        fxtemp = (yIC[i]*ri[i])/di[i]
        fx.append(fxtemp)
        
    i=0
    for i in range(num_bolts):     
        fytemp = (xIC[i]*ri[i])/di[i]
        fy.append(fytemp)
        
    i=0
    for i in range(num_bolts):    
        momenttemp = ri[i]*di[i]
        moment.append(momenttemp)
    
    totX = sum(fx)*0.75
    totY = sum(fy)*0.75
    totM = sum(moment)*0.75
    
    Pu_force = m.sqrt((totX*totX) + (totY*totY))
    
    Pu_moment = totM / (ecc - ICx*m.cos(m.radians(theta)) - ICy*m.sin(m.radians(theta)))
    
    res1 = [totX, totY, totM]
    
    res2 = [Pu_force, Pu_moment]
    
    return res1, res2
    
def ic_optimize(IC,xloc,yloc,ecc,theta):
    
    res1, res2 = ic_check(IC,xloc,yloc,ecc,theta)
    
    res = res2[0] - res2[1]
    
    return abs(res)
    
xloc = [-1.5,-1.5,-1.5,-1.5,1.5,1.5,1.5,1.5]
yloc = [-4.5,-1.5,1.5,4.5,4.5,1.5,-1.5,-4.5]

ICx = -1.69
ICy = 0.0
IC = [ICx,ICy]
theta = 0.0
ecc = 7.5
res1, res2 = ic_check(IC,xloc,yloc,ecc,theta)

test2 = ic_optimize(IC,xloc,yloc,ecc,theta)

test = sci.minimize(ic_optimize,IC,(xloc,yloc,ecc,theta), tol=1e-6)

IC_new = [test.x[0],test.x[1]]

res3, res4 = ic_check(IC_new,xloc,yloc,ecc,theta)

Open Source Structural Applications:
 
Upvote 0 Downvote
I improved the results slightly by also back checking sum Fx and sum Fy in my optimization function:
It is still very sensitive to the initial IC coordinate guess.

There is a great article on AISC for this by G. Donald Brandt from 2nd quarter 1982:
Python: 
def ic_optimize(IC,xloc,yloc,ecc,theta):
    
    res1, res2, check = ic_check(IC,xloc,yloc,ecc,theta)
    
    res = res2[0] - res2[1]
    
    px = res2[0]*m.sin(m.radians(theta))
    py = res2[0]*m.cos(m.radians(theta))
    
    fx = px + res1[0]
    fy = py - res1[1]
    
    res_magnitude = m.sqrt(res*res+fx*fx+fy*fy)
    
    return res_magnitude

Open Source Structural Applications:
 
Upvote 0 Downvote
Brandt's Method is much more reliable:

Python: 
from __future__ import division
import math as m
import scipy.optimize as sci

def ic_check(IC,xloc,yloc, ecc, theta):
    num_bolts = len(xloc)
    
    deltamax = 0.34
    Rult = 32.471
    
    ICx = IC[0]
    ICy = IC[1]
    xIC = []
    yIC = []
    di = []
    deltai = []
    ri = []
    
    fx = []
    fy = []
    moment = []
    
    for x in xloc:
        xICtemp = x - ICx
        xIC.append(xICtemp)
        
    for y in yloc:
        yICtemp = y - ICy
        yIC.append(yICtemp)
    
    i=0
    for i in range(num_bolts):
        ditemp = m.sqrt((xIC[i]*xIC[i])+(yIC[i]*yIC[i]))
        di.append(ditemp)
    
    dmax = max(di)
    
    i=0
    for i in range(num_bolts):
        deltaitemp = (di[i]/dmax)*deltamax
        deltai.append(deltaitemp)
        
    i=0
    for i in range(num_bolts):
        ritemp = Rult*m.pow(1-m.pow(m.e,-10.0*deltai[i]),0.55)
        ri.append(ritemp)
        
    i=0
    for i in range(num_bolts):    
        fxtemp = (yIC[i]*ri[i])/di[i]
        fx.append(fxtemp)
        
    i=0
    for i in range(num_bolts):     
        fytemp = (xIC[i]*ri[i])/di[i]
        fy.append(fytemp)
        
    i=0
    for i in range(num_bolts):    
        momenttemp = ri[i]*di[i]
        moment.append(momenttemp)
    
    totX = sum(fx)*0.75
    totY = sum(fy)*0.75
    totM = sum(moment)*0.75
    
    Pu_force = m.sqrt((totX*totX) + (totY*totY))
    
    Pu_moment = totM / (ecc - ICx*m.cos(m.radians(theta)) - ICy*m.sin(m.radians(theta)))
    
    res1 = [totX, totY, totM]
    
    res2 = [Pu_force, Pu_moment]
    
    check_output = [xIC,yIC,di,deltai,ri,fx,fy]
    
    return res1, res2, check_output
    
def ic_optimize(IC,xloc,yloc,ecc,theta):
    
    res1, res2, check = ic_check(IC,xloc,yloc,ecc,theta)
    
    res = res2[0] - res2[1]
    
    px = res2[0]*m.sin(m.radians(theta))
    py = res2[0]*m.cos(m.radians(theta))
    
    fx = px - res1[0]
    fy = py - res1[1]
    
    res_magnitude = m.sqrt(res*res+fx*fx+fy*fy)
    
    return res_magnitude

def ic_brandt(IC, xloc, yloc, Mp):
    num_bolts = len(xloc)
    
    deltamax = 0.34
 
    ICx = IC[0]
    ICy = IC[1]
    xIC = []
    yIC = []
    di = []
    deltai = []
    ri = []
    
    fx = []
    fy = []
    moment = []
    
    for x in xloc:
        xICtemp = x - ICx
        xIC.append(xICtemp)
        
    for y in yloc:
        yICtemp = y - ICy
        yIC.append(yICtemp)
    
    i=0
    for i in range(num_bolts):
        ditemp = m.sqrt((xIC[i]*xIC[i])+(yIC[i]*yIC[i]))
        di.append(ditemp)
    
    dmax = max(di)
    
    i=0
    for i in range(num_bolts):
        deltaitemp = (di[i]/dmax)*deltamax
        deltai.append(deltaitemp)
        
    i=0
    for i in range(num_bolts):
        ritemp = m.pow(1-m.pow(m.e,-10.0*deltai[i]),0.55)
        ri.append(ritemp)
        
    i=0
    for i in range(num_bolts):    
        fxtemp = -1*(yIC[i]*ri[i])/di[i]
        fx.append(fxtemp)
        
    i=0
    for i in range(num_bolts):     
        fytemp = (xIC[i]*ri[i])/di[i]
        fy.append(fytemp)
        
    i=0
    for i in range(num_bolts):    
        momenttemp = ri[i]*di[i]
        moment.append(momenttemp)
    
    Mi = sum(moment)
    Rult = -1*Mp/Mi
    
    Rx = sum(fx) * Rult
    Ry = sum(fy) * Rult
    
    
    table = [xIC, yIC, di, deltai, ri, moment, fx, fy]
    
    return Rx, Ry, Mi, table
    
def brandt(xloc, yloc, P_xloc, P_yloc, P_angle):
    # Bolt Group Instantaneous Center using method by G. Donald Brandt
    # Rapid Determiniation of Ultimate Strength Of Eccentrically Loaded Bolt Groups
    # AISC Journal 1982 2nd Quarter
    
    detailed_output = []
    
    num_bolts = len(xloc)
    
    n = num_bolts
    
    detailed_output.append(num_bolts)
    
    #Bolt Group Centroid
    if len(xloc)<3:
        anchor_x_bar = (xloc[0]+xloc[1])/2.00
        anchor_y_bar = (yloc[0]+yloc[1])/2.00
    
    else:
        j=0
        x_tot=0
        y_tot=0
        
        for i in xloc:
            x_tot = x_tot+xloc[j]
            y_tot = y_tot+yloc[j]    
            j+=1
        
        anchor_x_bar = x_tot/len(xloc)
        anchor_y_bar = y_tot/len(yloc)
     
    cg_anchors = [anchor_x_bar, anchor_y_bar]
    detailed_output.append(cg_anchors)
    
    # J - Polar Moment of Inertial of Bolt Group
    # sum(x^2+y^2)
    sum_x_square = 0
    sum_y_square = 0
    
    i=0
    for i in range(num_bolts):
        sum_x_square = sum_x_square + (xloc[i]-anchor_x_bar)**2
        sum_y_square = sum_y_square + (yloc[i]-anchor_y_bar)**2
    
    J = sum_x_square + sum_y_square
    detailed_output.append(['J',J])
    
    Px = -1*m.cos(m.radians(P_angle))
    Py = -1*m.sin(m.radians(P_angle))
    
    detailed_output.append([Px,Py])
    
    Mo = (-1*Px*(P_yloc-anchor_y_bar))+(Py*(P_xloc-anchor_x_bar))
    
    detailed_output.append(Mo)
    
    ax = (-1*Py*J) / (n * Mo)
    ay = (Px*J) / (n*Mo)
    
    detailed_output.append([ax,ay])
    
    Mp = (-1*Px*(P_yloc-anchor_y_bar-ay))+(Py*(P_xloc-anchor_x_bar-ax))
    
    detailed_output.append(Mp)
    
    IC_initial = [anchor_x_bar+ax,anchor_y_bar+ay]
    
    Rx, Ry, Mi, table = ic_brandt(IC_initial,xloc,yloc, Mp)
    
    detailed_output.append([Rx, Ry, Mi, table,"First IC pass"])
    
    fxx = Px + Rx
    fyy = Py + Ry
    F = m.sqrt(fxx*fxx+fyy*fyy)
    
    detailed_output.append([fxx,fyy,F,"F"])
    
    ax_new = (-1*fyy*J)/(n*Mo)
    ay_new = (fxx*J) / (n*Mo)
    
    detailed_output.append(["ax",ax_new,"ay",ay_new])
    
    IC_new = IC_initial  
    
    count = 0
    iterations = 0
    while count<1000:
        
        IC_new = [IC_new[0]+ax_new,IC_new[1]+ay_new]
        Mp_new = (-1*Px*(P_yloc-IC_new[1]))+(Py*(P_xloc-IC_new[0]))
        
        Rx, Ry, Mi, table = ic_brandt(IC_new,xloc,yloc, Mp_new)
        
        fxx = Px + Rx
        fyy = Py + Ry
        F = m.sqrt(fxx*fxx+fyy*fyy)
        
        ax_new = (-1*fyy*J)/(n*Mo)
        ay_new = (fxx*J) / (n*Mo)
        
        if F <= 0.00001:
            iterations = count
            count = 1000          
            solution = 'yes'
        else:        
            count +=1
            solution = 'no'
    
    detailed_output.append([fxx,fyy,F])        
    detailed_output.append(IC_new)
    detailed_output.append([solution,iterations,count])
    
    detailed_output.append([Rx, Ry, Mi, table])
    
    
    
    Cu = abs(Mi/Mp_new)
    
    detailed_output.append([Mi,Mp_new,Cu])
    
    return detailed_output, IC_new, Cu

# From Scratch Method Testing Zone    
xloc = [-1.5,-1.5,-1.5,-1.5,1.5,1.5,1.5,1.5]
yloc = [-4.5,-1.5,1.5,4.5,4.5,1.5,-1.5,-4.5]

#Center of Anchor Group
if len(xloc)<3:
    anchor_x_bar = (xloc[0]+xloc[1])/2.00
    anchor_y_bar = (yloc[0]+yloc[1])/2.00

else:
    j=0
    x_tot=0
    y_tot=0
    
    for i in xloc:
        x_tot = x_tot+xloc[j]
        y_tot = y_tot+yloc[j]    
        j+=1
    
    anchor_x_bar = x_tot/len(xloc)
    anchor_y_bar = y_tot/len(yloc)

ICx = anchor_x_bar
ICy = anchor_y_bar
IC = [ICx,ICy]
theta = 0
ecc = 4
res1, res2, check1 = ic_check(IC,xloc,yloc,ecc,theta)

test2 = ic_optimize(IC,xloc,yloc,ecc,theta)

test = sci.minimize(ic_optimize,IC,(xloc,yloc,ecc,theta),method='Nelder-Mead', tol=1e-6)

IC_new = [test.x[0],test.x[1]]

res3, res4, check2 = ic_check(IC_new,xloc,yloc,ecc,theta)
    
C = res4[0]/(0.75*32.471)   

# Brandt's Method Testing Zone
# angle is positive counter-clockwise and measured from 0 pointing to the right
# anchor and P coordinates are global
x_b = [-1.5,-1.5,-1.5,-1.5,1.5,1.5,1.5,1.5]
y_b = [-4.5,-1.5,1.5,4.5,4.5,1.5,-1.5,-4.5]
P_xloc = 7
P_yloc = 0
P_angle = 90

brandt = brandt(x_b, y_b, P_xloc, P_yloc, P_angle)

Open Source Structural Applications:
 
Upvote 0 Downvote
thank you a ton Celt.

A few questions for you as I am trying to parse this code.

What does return do in python:
Python: 
        momenttemp = ri[i]*di[i]
        moment.append(momenttemp)
    
    totX = sum(fx)*0.75
    totY = sum(fy)*0.75
    totM = sum(moment)*0.75
    
    Pu_force = m.sqrt((totX*totX) + (totY*totY))
    
    Pu_moment = totM / (ecc - ICx*m.cos(m.radians(theta)) - ICy*m.sin(m.radians(theta)))
    
    res1 = [totX, totY, totM]
    
    res2 = [Pu_force, Pu_moment]
    
    check_output = [xIC,yIC,di,deltai,ri,fx,fy]
    
    [b][u]return[/u][/b] res1, res2, check_output

It looks like the program now knows the variables res1, res2, and check_output after the function is performed. Is this correct?

Could the (2) functions you have created be combined into one function? It looks like ic_optomize just builds on results from ic_check? The way I am interrupting this, res_magnitude (the final output from ic_optomize) could be programmed into ic_check.

Could the minimize function then be run ic_check in this case? Or does the function that is being minimized only allowed to have (1) output?


S&T
 
Upvote 0 Downvote
i will also dig into the Brandt method when I get back to the office to login into AISC.

Looks like nice method for quick convergence.

S&T
 
Upvote 0 Downvote
The return line at the end of the function specifies what variable or value can be assigned as the function result to be used globally for further operations.

I have limited knowledge of how to properly use the scipy.optimize.minimize function but it seemed to want one result to check against, i also found that if the result was negative it would keep going negative rather than adjust to start returning back closer to 0.

Open Source Structural Applications:
 
Upvote 0 Downvote
Cool stuff Celt! Digging into the GUI stuff this weekend once I get some free time.

S&T
 
Upvote 0 Downvote
I have now set this up to run using the Scipy.optimize root function:

Python: 
def xlic_calc(IC, xloc, yloc, vals):
    ecc =vals[0]
    theta = np.radians(vals[1])
    deltamax = vals[2]
    Rult = vals[3]
    num_bolts = len(xloc)
    
    ICx = IC[0]
    ICy = IC[1]
    Pu = IC[2]
    Pux = Pu * np.sin(theta)
    Puy = Pu * np.cos(theta)
    Mu = Pu * (ecc-ICx*np.cos(theta)-ICy*np.sin(theta))
    xIC = []
    yIC = []
    di = []
    deltai = []
    ri = []
    
    fx = []
    fy = []
    moment = []
    
    for x in xloc:
        xICtemp = x - ICx
        xIC.append(xICtemp)
        
    for y in yloc:
        yICtemp = y - ICy
        yIC.append(yICtemp)
    
    i=0
    for i in range(num_bolts):
        ditemp = np.sqrt((xIC[i]*xIC[i])+(yIC[i]*yIC[i]))
        di.append(ditemp)
    
    dmax = max(di)
    
    i=0
    for i in range(num_bolts):
        deltaitemp = (di[i]/dmax)*deltamax
        deltai.append(deltaitemp)
    
        ritemp = Rult*(1-np.e**(-10.0*deltai[i]))**0.55
        ri.append(ritemp)
                    
        fxtemp = (yIC[i]*ri[i])/di[i]
        fx.append(fxtemp)
             
        fytemp = (xIC[i]*ri[i])/di[i]
        fy.append(fytemp)
          
        momenttemp = ri[i]*di[i]
        moment.append(momenttemp)
    
    totX = sum(fx)*0.75
    totY = sum(fy)*0.75
    totM = sum(moment)*0.75
    
    return [totX, totY, totM, Pux, Puy, Mu]

It's set up to be called from Excel using xlwings, but you can also call direct from Python:

>>> res = sopt.root(ic_check, x0, vals2)
>>> res.x
array([-1.68496865e+00, 4.01869965e-16, 7.53006559e+01])

Doug Jenkins
Interactive Design Services
 
Upvote 0 Downvote
continued:
Inputs are:
x0, the initial guess: [-3, 0,100]
vals2: [[target vals],[[Xloc],[Yloc],[knowns]]]

[target vals] = [0,0,0]
[knowns] = [Eccentricity, Theta, deltamax,Rult]

I didn't have any problem with sensitivity to the initial guess values, and it works very fast with the default solver. Several of the other solvers work very slowly, or not at all.

I will post the Excel version in the next few days. Here's what it looks like:

Scipy-solve_l9l10n.png



Doug Jenkins
Interactive Design Services
 
Upvote 0 Downvote
Part 3.

The ic_check function, which is the one you call from sopt.root:

Python: 
def ic_check(IC, vals):   
    xloc = vals[1][0]
    yloc  = vals[1][1]
    vals2 =vals[1][2]
    
    res1 = xlic_calc(IC,xloc, yloc, vals2)
    
    totX = res1[0]
    totY = res1[1]
    totM = res1[2]
            
    Pu =IC[2]
    Pux = res1[3]
    Puy = res1[4]
    Mu = res1[5]
    
    res3 = [totX-Pux, totY-Puy, totM-Mu]
    
    return res3

Doug Jenkins
Interactive Design Services
 
Upvote 0 Downvote
@IDS:
I run into some problems with loads at flatter angles and small eccentricities, do you get the same issue. As an example try the load at an eccentricity of 5 and angle of 15 with that same 8 bolt layout.

Capture_vmpojo.png

Capture1_yma6ug.png


Open Source Structural Applications:
 
Upvote 0 Downvote
It was a tolerance issue certain bolt and load configurations don't converge beyond a certain precision. I added the ability to control the tolerance.

the 8 bolt group needed a tolerance of 0.00001
a 3 bolt group (0,0) (0,3) (0,6) load at (5,3) angle: 15 needed a tolerance of 0.1 on my end

I'm sure the scipy module may have something built in to handle these odd conditions.

Open Source Structural Applications:
 
Upvote 0 Downvote
Status
Not open for further replies.

Similar threads

MechEngineer2012
  • Question Question
Python programming assistance needed
Replies
2
Views
5K
IRstuff
IRstuff
GregLocock
  • Locked
  • Question Question
New integrated AI/Python browser based app - Dystr
Replies
3
Views
4K
Sparweb
Sparweb
ali-gct-meca
  • Locked
  • Question Question
Hello, I have a pressure vessel
Replies
0
Views
759
ali-gct-meca
ali-gct-meca
KA786
  • Locked
  • Question Question
I need help in sympy programming in Python any expert advise most welcome !!!!
Replies
2
Views
381
IRstuff
IRstuff
Onda
  • Locked
  • Question Question
MSC.Software Apex - Dedicated Forum
Replies
4
Views
4K
Onda
Onda

Part and Inventory Search

Sponsor

Back
Top