Tugas Besar Metode Numerik - Josiah Enrico S

From ccitonlinewiki
Revision as of 13:13, 14 January 2021 by JosiahEnrico (talk | contribs) (Result and Analysis)
Jump to: navigation, search

Abstract

Introduction

Objective

  • Membuat program analisis truss dan optimasi sederhana dengan Modelica
  • Mengoptimasi harga pembuatan rangka truss sederhana dengan memvariasi dimensi dan elastisitas material.

Methodology

  • Metode Optimasi Golden Section (univariabel)
  • Metode Optimasi Golden Section (multivariabel)

Procedure

Geometri dan Load


Force
F1 2000
F2 1000

Constraint:

- Spesifikasi L (Panjang) dan geometri rangka truss

- Gaya beban terhadap struktur (1000 N dan 2000 N)


Asumsi:

- Beban akan terdistribusi hanya pada point penghubung dan semua gaya yang dipengaruhi moment diamggap tidak ada karena sistem bersifat truss


Koleksi Data

Data Proceessing

Truss Structure Analysis Model

Pada tahap pertama, dilakukanka simulasi menggunakan Modelica untuk mengetahui pengaruh beban gaya yang diberikan terhadap struktur truss. Pengaruh ini direpresantasikan oleh 4 nilai yaitu displacement/perpindahan, gaya reaksi, tekanan di setiap truss dan faktor keamanan di setiap truss. Keempat nilai ini diperoleh dengan menggunakan operasi tranformasi matriks untuk pembebanan truss dalam Finite Element Analysis. Operasi ini diterjemahkan ke dalam bahasa C dalam pemrograman Modelica. Berikut ini adalah hasil simulasi Modelica dan program tertulisnya:

model Truss_3D_Tugas_Besar_Safety
// Base Code by Josiah Enrico S. 1906356286
// Editing and Testing by Ahmad M. Fahmi 1806181836
// Commenting and Cleaning by Christopher S. Erwin 1806201056 

// === DEFINITION PHASE ===

// Define the properties of the truss arrangement
parameter Integer Nodes = size(node_points,1);    // Number of Nodes
parameter Integer Members = size(node_pairs,1);   // Number of Members

// Define the properties of the truss material
parameter Real Yield = 215e6;           // Yield Strength (Pa)
parameter Real Area = 0.000224;         // Area L Profile [Dimension=0.03, Thickness=0,004] (m^2)
parameter Real Elas = 193e9;            // Elasticity of material [SS 304] (Pa)

// Define connections; list pair of nodes joined by a member
parameter Integer node_pairs[:,2]=[1,5; 2,6; 3,7; 4,8;         // Lowest vertical members
                                  5,6; 6,7; 7,8; 5,8;         // 1st shelf
                                  5,9; 6,10; 7,11; 8,12;      // Middle vertical members
                                  9,10; 10,11; 11,12; 9,12;   // 2nd shelf
                                  9,13; 10,14; 11,15; 12,16;  // Top vertical members
                                  13,14; 14,15; 15,16; 13,16];// 3rd shelf
                                                            
// Define coordinates and boundaries of each node
// [x, y, z, boundary_x, boundary_y, boundary_z]
// If the node is bound (stationary) along that axis, then boundary_axis = 1
// otherwise the node is free and boundary_axis = 0
parameter Real node_points[:,6]=[0.3,-0.375,0,1,1,1;      // Node 1
                                -0.3,-0.375,0,1,1,1;     // Node 2
                                -0.3,0.375,0,1,1,1;      // Node 3
                                0.3,0.375,0,1,1,1;       // Node 4
                      
                                0.3,-0.375,0.6,0,0,0;    // Node 5
                                -0.3,-0.375,0.6,0,0,0;   // Node 6
                                -0.3,0.375,0.6,0,0,0;    // Node 7
                                0.3,0.375,0.6,0,0,0;     // Node 8
                           
                                0.3,-0.375,1.2,0,0,0;    // Node 9
                                -0.3,-0.375,1.2,0,0,0;   // Node 10  
                                -0.3,0.375,1.2,0,0,0;    // Node 11
                                0.3,0.375,1.2,0,0,0;     // Node 12
                      
                                0.3,-0.375,1.8,0,0,0;    // Node 13
                                -0.3,-0.375,1.8,0,0,0;   // Node 14
                                -0.3,0.375,1.8,0,0,0;    // Node 15
                                0.3,0.375,1.8,0,0,0];    // Node 16
                           
// Define external forces in Newtons on each node
// [Fx, Fy, Fz]
parameter Real Forces[nodes_xyz]={0,0,0,      // Node 1
                                 0,0,0,      // Node 2
                                 0,0,0,      // Node 3
                                 0,0,0,      // Node 4
                                 0,0,0,      // Node 5
                                 0,0,0,      // Node 6
                                 0,0,0,      // Node 7
                                 0,0,0,      // Node 8
                                 0,0,0,      // Node 9
                                 0,0,0,      // Node 10
                                 0,0,0,      // Node 11
                                 0,0,0,      // Node 12
                                 0,0,-500,   // Node 13
                                 0,0,-1000,  // Node 14
                                 0,0,-1000,  // Node 15
                                 0,0,-500};  // Node 16

// === SOLUTION PHASE ===

Real displacement[nodes_xyz], reaction[nodes_xyz];
Real check[3]; // To check if forces in xyz directions sum to zero (structure is static)

Real safety[Members];     // Safety Factor of each member (based on stress experienced)
Real mem_dis[3];          // Displacement in xyz directions experienced by member
Real stress_vector[3];    // Stress experienced by member in xyz directions
Real stress_mag[Members]; // Magnitude of stress experienced by member

protected
parameter Integer nodes_xyz = 3 * Nodes;  // Size of vector of {node_#_x, node_#_y, node_#_z} containing all nodes
Real node_i[3], node_j[3];                // Coordinates in xyz for member's node pairs 

// Global Stiffness Matrix construction
Real G_element[nodes_xyz, nodes_xyz], G_total[nodes_xyz, nodes_xyz], G_total_new[nodes_xyz, nodes_xyz];

// Transformation Matrix construction
Real cos_x, cos_y, cos_z, L, T[3,3];

// Limits to eliminate Float Errors
Real err_solve = 10e-10, err_check = 10e-4;

algorithm
//Creating Global Matrix
G_total := identity(nodes_xyz);

for mem in 1:Members loop
 // Get coordinates of member's nodes
 for i in 1:3 loop
   node_i[i] := node_points[node_pairs[mem,1],i];
   node_j[i] := node_points[node_pairs[mem,2],i];
 end for;
 
 // Finding the Stiffness Matrix [Ke]
 // {F} = [Ke]{U}
 // [Ke] = k[T]
 L := Modelica.Math.Vectors.length(node_j - node_i);
 cos_x := (node_j[1]-node_i[1])/L;
 cos_y := (node_j[2]-node_i[2])/L;
 cos_z := (node_j[3]-node_i[3])/L; 
 T := (Area * Elas / L)*[cos_x^2,cos_x*cos_y,cos_x*cos_z;
                         cos_y*cos_x,cos_y^2,cos_y*cos_z;
                         cos_z*cos_x,cos_z*cos_y,cos_z^2];

 // Transforming to Global Matrix
 G_element := zeros(nodes_xyz,nodes_xyz); 
 for m,n in 1:3 loop
   G_element[3*(node_pairs[mem,1]-1)+m, 3*(node_pairs[mem,1]-1)+n] :=  T[m,n];
   G_element[3*(node_pairs[mem,2]-1)+m, 3*(node_pairs[mem,2]-1)+n] :=  T[m,n];
   G_element[3*(node_pairs[mem,2]-1)+m, 3*(node_pairs[mem,1]-1)+n] := -T[m,n];
   G_element[3*(node_pairs[mem,1]-1)+m, 3*(node_pairs[mem,2]-1)+n] := -T[m,n];
 end for;
 
 // Adding to total Global Matrix
 G_total_new := G_total + G_element;
 G_total := G_total_new;
end for;

//Implementing boundary conditions to Global Matrix
for nod in 1:Nodes loop
 // Boundary x-axis
 if node_points[nod,4] <> 0 then
   for i in 1:Nodes*3 loop
     G_total[(nod*3)-2,i] := 0;
     G_total[(nod*3)-2,(nod*3)-2] := 1;
   end for;
 end if;
 
 // Boundary y-axis
 if node_points[nod,5] <> 0 then
   for i in 1:Nodes*3 loop
     G_total[(nod*3)-1,i] := 0;
     G_total[(nod*3)-1,(nod*3)-1] := 1;
   end for;
 end if;
 
 // Boundary z-axis
 if node_points[nod,6] <> 0 then
   for i in 1:Nodes*3 loop
     G_total[nod*3,i] := 0;
     G_total[nod*3,nod*3] := 1;
   end for;
 end if;
end for;

// Solving displacement (U)
// [G_total]{U} = {F}
displacement := Modelica.Math.Matrices.solve(G_total, Forces);

// Solving reaction (R)
// {R} = [G_total]{U} - {F}
reaction := (G_total_new * displacement) - Forces;

// Eliminating float errors
for i in 1:nodes_xyz loop
reaction[i] := if abs(reaction[i])<= err_solve then 0 else reaction[i];
displacement[i] := if abs(displacement[i])<= err_solve then 0 else displacement[i];
end for;

// Checking Sum of Forces
// For static structures, should add up to 0
check[1]:=sum({reaction[i] for i in (1:3:(nodes_xyz-2))}) + sum({Forces[i] for i in (1:3:(nodes_xyz-2))});
check[2]:=sum({reaction[i] for i in (2:3:(nodes_xyz-1))}) + sum({Forces[i] for i in (2:3:(nodes_xyz-1))});
check[3]:=sum({reaction[i] for i in (3:3:(nodes_xyz))}) + sum({Forces[i] for i in (3:3:(nodes_xyz))});

// Eliminating float errors
for i in 1:3 loop
 check[i] := if abs(check[i])<= err_check then 0 else check[i];
end for;

// === POSTPROCESSING PHASE ===

//Calculating stress in each truss
for mem in 1:Members loop
 // Getting node coordinate data and displacement for current member
 for i in 1:3 loop
   node_i[i] := node_points[node_pairs[mem,1],i];
   node_j[i] := node_points[node_pairs[mem,2],i];
   mem_dis[i] := abs(displacement[3*(node_pairs[mem,1]-1)+i] - displacement[3*(node_pairs[mem,2]-1)+i]);
 end for;
 
 // Solving the Displacement Matrix for current member
 // Stress = F/A = (E/L) * [dL]
 L:=Modelica.Math.Vectors.length(node_j-node_i);
 cos_x:=(node_j[1]- node_i[1])/L;
 cos_y:=(node_j[2]- node_i[2])/L;
 cos_z:=(node_j[3]- node_i[3])/L; 
 T := (Elas/L)*[cos_x^2,cos_x*cos_y,cos_x*cos_z;
               cos_y*cos_x,cos_y^2,cos_y*cos_z;
               cos_z*cos_x,cos_z*cos_y,cos_z^2];
 
 stress_vector:=(T*mem_dis); // The stress in vector components {stress_x, stress_y, stress_z}
 stress_mag[mem]:=Modelica.Math.Vectors.length(stress_vector); // Magnitude/length of the stress vectors for each member
end for;

// Safety Factor
// Factor of Safety = (Yield Strength of material) / (Stress experienced by Member)
for mem in 1:Members loop
 if stress_mag[mem]>0 then
   safety[mem] := Yield / stress_mag[mem];
 else
   safety[mem] := 0;
 end if; 
end for;

end Truss_3D_Tugas_Besar_Safety;


Cost and Cost Efficiency of Area Variation and Constrained Elasticity

Berdasarkan hasil analisis dari struktur truss sebelumnya, pada bagian ini parameter area cross-section dari setiap truss dalam struktur tesebut divariasi dengan menggunakan yann standar ukuran yang tertera pada prosedur. Data perbedaan Area dalam material yang sama tersebut diproses menjadi variasi harga (cost) dan efisiensi harga (safaty factor/cost) agar dapat diproses di akhir dengan optimisasi metode Powell. Dalam Prosesnya, data harga material yang terbatas di pasar akan di-regresi orde 2 untuk mendapatkan nilai harga pada ukuran-ukuran yang tidak memilik harga. Berikut adalah tabel hasil proses data dan grafik yang dihasilkan:

Cost and Cost Efficiency of Elasticity Variation and Constrained Area

Mirip dengan pengolahan data sebelumnya, bagian ini juga akan menghasilkan nilai variasi harga (cost) dan efisiensi harga (safaty factor/cost) untuk diolah dengan optimisasi metode powell di akhir. Namun, berbeda dengan sebalumnya, pada bagian ini area cross section akan disamakan nilainya dan elastisitas material akan divariasikan. Dalam pengolah ini, metode curve fitting orde 2 juga dilakukan untuk mendapatkan nilai-nilai harga, kekuatan yield, densitas yang kosong. Berikut adalah tabel hasil proses data dan grafik yang dihasilkan:

Result and Analysis

Hasil dari studi ini adalah nilai cost efficincy optimal yang dipengaruhi oleh 2 faktor yaitu elastisitas dan area. Oleh karena itu jenis optimisasi yang dilakukakn adalah optimasi multivariabel (lebih dari satu variabel) dengan metode powell. Metode ini dipilih karena tergolong cukup sederhana dan cepat mencapai titik konvergen. Secara singkat, agar lebih memahami optimasi multivariabel ini, bisa dibayangkan distribusi cost efficieny digambarkan sebagi grafik 3 dimensi dalam sumbu x, y, z, Cost efficiency akan diwakili skala z, elastisitas sebagai skala y dan area sebagai skala x. Nilai dari cost efficincy(z) pada suatu titik adalah fungsi area(x) dan elastisitas (y). Persamaan kuadrat cost efficiency yang diperolah dari curve fitting tahap sebelumnya pada 2 variasi, yaitu area dan elastisitas, akan di rata-ratakan sehingga menjadi fungsi tunggal. Hasil tersebut akan menghasikan permukaan kontur bidang xy yang memiliki titik optimum. koordinat titik tersebut akan dicari secara linear satu persatu dengan metode Golden. Hasil dari simulasi dan program tertulis Modelica akan ditampilkan sebagai berikut:

Modelica Result

Powell Method

model Powell_Method_Opt_Jos_Ratio

parameter Integer order=2; //order of regression
parameter Integer GoldN=15; //maximum iteration of linear gold optimization
parameter Integer N=10; //maximum iteration
parameter Real maxerror=1e-50; //maximum error

//Assumed to be X function
parameter Real RatioArea[size(Area,1)]={3.6534e-5, 3.75348e-5, 3.79639e-5, 3.79507e-05, 3.74821e-05};
parameter Real CostArea[size(Area,1)]={829776,1053955,1308215,1722235,2202651};
parameter Real Area[:]={141e-6, 184e-6, 231e-6, 304e-6, 384e-6};
Real CoeArea[order+1];
Real CoeAreaCost[order+1];

//Assumed to be Y function
parameter Real RatioElas[size(Elas,1)]={2.83875e-5, 8.5637e-5, 12.5153e-5, 5.70447e-5};
parameter Real CostElas[size(Elas,1)]={1367994,732878,492691,622939};

parameter Real Elas[:]={192e9, 195e9, 198e9, 201e+9};
Real CoeElas[order+1];
Real CoeElasCost[order+1];

//Guessing Start (1 point, 2 vector) - the second vector will be auto-generated
parameter Real StartPoint[2]={171e-6,195e9};
parameter Real VectorPoint[2]={205e-6,200e9};

//Solution
Real RatioZ[N]; //Optimum result
Real CostZ; //Cost Optimum result
Real x3[2]; //Optimum point coordinate (x,y)
Real error[N]; //Error 

protected
Real jos, RatioZX, RatioZY, CostZX, CostZY;
Real x0tes[2], x0[2], x1[2], x2[2];
Real grad1, grad2, grad3;

//Checking approach 
Real cek0[N], cek1[N], cek2[N], cek3[N], yek0[N], yek1[N], yek2[N], yek3[N];

algorithm
//Creating function at both 
CoeArea:=Curve_Fitting(Area,RatioArea,order);
CoeElas:=Curve_Fitting(Elas,RatioElas,order);
CoeAreaCost:=Curve_Fitting(Area,CostArea,order);
CoeElasCost:=Curve_Fitting(Elas,CostElas,order);

//Creating 2 vector defined in 2 gradient (grad1,grad2)
x0:=StartPoint;
x0tes:=VectorPoint;
grad1:=(x0tes[2]-x0[2])/(x0tes[1]-x0[1]);
grad2:=-0.5*grad1;

//Provide error variable
RatioZX:=CoeArea[1]*x0[1]^2+CoeArea[2]*x0[1]+CoeArea[3];
RatioZY:=CoeElas[1]*x0[2]^2+CoeElas[2]*x0[2]+CoeElas[3];
jos:=(RatioZX+RatioZY)/2;

//Computing approach
for i in 1:N loop
 cek0[i]:=x0[1];
 yek0[i]:=x0[2];
 x1:=Gold_Opt_Func_2D_Vector(x0,grad1,GoldN,CoeArea,CoeElas);
 cek1[i]:=x1[1];
 yek1[i]:=x1[2];
 x2:=Gold_Opt_Func_2D_Vector(x1,grad2,GoldN,CoeArea,CoeElas);
 cek2[i]:=x2[1];
 yek2[i]:=x2[2];
 
   if (x0[1]==x2[1]) and (x0[2]<>x2[2]) then
     grad3:=1/0; //divergent gradient, please subtitute vector
     x3:=Gold_Opt_Func_2D_Vector_1(x2,1,GoldN,CoeArea,CoeElas);    
   else  
     grad3:=(x2[2]-x0[2])/(x2[1]-x0[1]);
     x3:=Gold_Opt_Func_2D_Vector(x0,grad3,GoldN,CoeArea,CoeElas);  
   end if;
 
 cek3[i]:=x3[1];
 yek3[i]:=x3[2];
 grad1:=grad2;
 grad2:=grad3; 
   
 //Function Created
 RatioZX:=CoeArea[1]*x3[1]^2+CoeArea[2]*x3[1]+CoeArea[3];
 RatioZY:=CoeElas[1]*x3[2]^2+CoeElas[2]*x3[2]+CoeElas[3];
 RatioZ[i]:=(RatioZX+RatioZY)/2;
 
 CostZX:=CoeAreaCost[1]*x3[1]^2+CoeAreaCost[2]*x3[1]+CoeAreaCost[3];
 CostZY:=CoeElasCost[1]*x3[2]^2+CoeElasCost[2]*x3[2]+CoeElasCost[3];
 CostZ:=(CostZX+CostZY)/2;
 x0:=x3;
 
 error[i]:=abs(1-jos/RatioZ[i]);
 jos:=RatioZ[i];
 if error[i]<maxerror then
   break;
 end if;
end for;

end Powell_Method_Opt_Jos_Ratio;

Gold Optimization Function

function Gold_Opt_Func_2D_Vector

//only for second order polynomial
input Real Point[2];
input Real gradient; //Gradient
input Integer N; //Maximum iteration
input Real coe1[3];
input Real coe2[3];
output Real optimum[2];

protected
Real C=Point[2]-Point[1]*gradient;
Real xlo=141e-6;
Real xhi=384e-6;

Real f1[N], f2[N], x1[N], x2[N], ea[N];
Real fx;
Real es=1e-10; // maximum error
Real d, xl, xu, xint, R=(5^(1/2)-1)/2;

algorithm
xl := xlo; 
xu := xhi;
 
for i in 1:N loop
 d:= R*(xu-xl);
 x1[i]:=xl+d;
 x2[i]:=xu-d;

 f1[i]:=((coe1[1]*x1[i]^2+coe1[2]*x1[i]+coe1[3])+(coe2[1]*(gradient*x1[i]+C)^2+coe2[2]*(gradient*x1[i]+C)+coe2[3]))/2;
 f2[i]:=((coe1[1]*x2[i]^2+coe1[2]*x2[i]+coe1[3])+(coe2[1]*(gradient*x2[i]+C)^2+coe2[2]*(gradient*x2[i]+C)+coe2[3]))/2;
 xint:=xu-xl;
 
 if f1[i]>f2[i] then
   xl:=x2[i];
   optimum[1]:=x1[i];
   optimum[2]:=(gradient*x1[i]+C);
   fx:=f1[i];
   else
     xu:=x1[i];
     optimum[1]:=x2[i];
     optimum[2]:=(gradient*x2[i]+C);
     fx:=f2[i];
 end if;
 
 ea[i]:=(1-R)*abs((xint)/optimum[1]);
 if ea[i]<es then
   break;
 end if;
end for;

end Gold_Opt_Func_2D_Vector;

Curve Fitting

function Curve_Fitting

input Real X[:];
input Real Y[size(X,1)];
input Integer order=2;
output Real Coe[order+1];

protected
Real Z[size(X,1),order+1];
Real ZTr[order+1,size(X,1)];
Real A[order+1,order+1];
Real B[order+1];

algorithm

for i in 1:size(X,1) loop
 for j in 1:(order+1) loop
 Z[i,j]:=X[i]^(order+1-j);
 end for;
end for;
ZTr:=transpose(Z);

A:=ZTr*Z;
B:=ZTr*Y;
Coe:=Modelica.Math.Matrices.solve(A,B);

end Curve_Fitting;

Acknowledgement

Reference