Eng-Tips is the largest engineering community on the Internet

Intelligent Work Forums for Engineering Professionals

Register Log in
  • Congratulations waross 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
Apr 7, 2015
468
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
 
Replies continue below

Recommended for you

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
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

ali-gct-meca
  • Locked
  • Question
Hello, I have a pressure vessel
Replies
0
Views
36
ali-gct-meca
ali-gct-meca
KA786
  • Locked
  • Question
I need help in sympy programming in Python any expert advise most welcome !!!!
Replies
2
Views
21
IRstuff
IRstuff
CET developer
  • Locked
  • Question
Using CET.SteelConnDesign V1.5 and FreeCAD for Automated 2D Drawing Generation
Replies
0
Views
53
CET developer
CET developer
enghelp2806
  • Locked
  • Question
Anyone recognize this FEA solver language?
Replies
8
Views
40
FEMAP_MAN
FEMAP_MAN
peter_engiio
  • Locked
  • Question
Help building a companion app for Risa
Replies
0
Views
15
peter_engiio
peter_engiio

Part and Inventory Search

Sponsor

Back
Top