UAS

Perhitungan Optimasi Struktur

1. Truss Element

import numpy as np def meshtruss (p1,p2,nx,ny):

   nodes = []
   bars = []
   xx = np.linspace(p1[0],p2[0],nx+1)
   yy = np.linspace(p1[1],p2[1],ny+1)
   for y in yy:
       for x in xx :
           nodes.append([x,y])

   for j in range (ny):
       for i in range (nx):
           n1 = i+j*(nx+1)
           n2 = n1+1
           n3 = n1+nx+1
           n4 = n3 + 1
           bars.extend([[n1,n2],[n1,n3],[n1.n4],[n2,n3]])
       bars.append ([n2,n4])
   index = ny*(nx+1)+1
   for j in range(nx):
       bars.append([index+j-1,index+j])
   return np.array(nodes), np.array(bars)

from openopt import NLP from matplotlib.pyplot import figure,show

def opttruss (coord,connec,E,F,freenote,V0,plotdisp=False):

   n= connec.shape[0]
   m= coord.shape[0]
   vector = coord[connec[:,1],:]-coord[connec[:,0],:]
   l = np.sqrt((vectors**2).sum(axis=1))
   e = vectors.T/1
   B = (e[np.newaxis]*e[:,np.newaxis]).T

   def fobj(x):
       D = E*x/1
       kx = e*D
       K = np.zeros((2*m,2*m))
       for i in range (n):
           aux = 2*connec[i,:]
           index = np.r_[aux[0]:aux[0]+2,aux[1]:aux[1]+2
           k0 = np.concatenate((np.concatenate((B[i],-B[i], axis=1),no.concatenate ((-B[i],B[i],axis=1)),axis=0)
           K[np.ix_(index,index)]=K[np.ix_(index,index)] + D[i]* k0

       block = freenote.flatten().nonzero()[0]
       matrix = K[np.ix_(block,block)]
       rhs = F.flatten()[block]
       solution = np.linalg.solve(matrix,rhs)
       u = freenode.astype('float').flatten()
       u[block] = solution

dst...

2. Golden Section Search

The golden section search is the counterpart of bisection used in finding roots of equations. Suppose that the minimum of f(x) has been bracketed in the interval (a,b) of length h. To telescope the interval, we evaluate the function at x1 =b− Rh and x2 =a + Rh, as shown in Figure 10.2(a). The constant R is to be determined shortly. If f1 > f2 as indicated in the figure, the minimum lies in (x1,b); otherwise itislocatedin(a, x2). Assuming that f1 > f2, we set a ←x1 and x1 ←x2, which yields a new interval (a,b) of length h� = Rh, as illustrated in Figure 10.2(b). To carry out the next telescopingoperationweevaluatethefunctionat x2 =a + Rh� andrepeattheprocess. The procedure works only if Figures 10.1(a) and (b) are similar (i.e., if the same constant R locatesx1 andx2 inbothfigures).ReferringtoFigure10.2(a),wenotethat x2 −x1 =2Rh−h.ThesamedistanceinFigure10.2(b)isx1 −a =h�− Rh�.Equating thetwo,we get

2Rh−h=h�− Rh�

code : For Profil L import math def bracket(f,x1,h):

   c = 1.618033989
   f1 = f(x1)
   x2 = x1 + h
   f2 = f(x2)
   if f2 > f1:
       h = -h
       x2 = x1 + h
       f2 = f(x2)
       if f2 > f1:
           return x2,x1 - h
   for i in range (100):
       h = c*h
       x3 = x2 + h
       f3 = f(x3)
       if f3 > f2:
           return x1,x3
       x1 = x2
       x2 = x3
       f1 = f2
       f2 = f3
       print ("Bracket did not find a minimum")

def search(f,a,b,tol=1.0e-9):

       nIter = int(math.ceil(-2.078087*math.log(tol/abs(b-a))))
       R = 0.618033989
       C = 1.0 - R
       x1 = R*a + C*b
       x2 = C*a + R*b
       f1 = f(x1)
       f2 = f(x2)
       for i in range(nIter):
           if f1 > f2:
               a = x1
               x1 = x2
               f1 = f2
               x2 = C*a + R*b
               f2 = f(x2)
           else:
               b = x2
               x2 = x1
               f2 = f1
               x1 = R*a + C*b
               f1 = f(x1)
       if f1 < f2:
           return x1,f1
       else:
           return x2,f2

print("Aplikasi Optimasi Section Modulus L Stiffner") print("Kondisi Terikat : lebar alas > lebar atas > lebar tengah") b1 = eval(input("Nilai lebar bangun alas (meter) :")) b3 = eval(input("Nilai lebar bangun atas (meter):")) b2 = eval(input("Nilai lebar bangun tengah (meter) :")) H = eval(input("Nilai tinggi T stiffner (meter) :")) def f(x):

   A1 = b1*(H-x)/2
   A2 = b2*x
   A3 = b3*(H-x)/2
   d1 = 1/2*(H-x)/2
   d2 = 1/2*x+(H-x)/2
   d3 = 3/4*(H-x)+x
   I1 = 1/12*b1*((H-x)/2)**3
   I2 = 1/12*b2*x**3
   I3 = 1/12*b3*((H-x)/2)**3
   dc = H-(d1*A1+d2*A2+d3*A3)/(A1+A2+A3)
   I = I1-A1*(d1-dc)**2+I2-A2*(d2-dc)**2+I3-A3*(d3-dc)**2
   Z = I/dc
   return Z

xStart = 0.0 h = 1.0 x1,x2 = bracket(f,xStart,h) y,fMin = search(f,x1,x2) print("optimal sectional area =",-fMin) print("sectional area awal" , f(H)) A = -fMin/f(H)*100 print ("efisiensi",A,"%") input ("\nPress return to exit")