I'm trying develop a structural weight optimization of a Glass Reinforcement Plastic (GRP) High Speed Craft using MINLP and the GEKKO package. I need to consider many design variables as integers and other as real and continuous variables. To evaluate the constraints, I'm following the requirements of ISO 12215 part 5. This rule establishes the minimum structural requirements for a craft. The problem is that I need use many conditional like (if, max or min), I read that GEKKO changes these types of conditional statements because the traditional conditions are not continuous. Therefore, I use m.if3, m.max3, or m.min3. However, when I run the program these GEKKO conditional statements do not work. Many operations (intermediates) take infinity value because the program divides by zero and never converges to a solution. Any suggestions to understand the solver and obtain a successful solution?
My problem is relative short because has 20 design variables with around 15 constraints without time dependence.
# Optimization (Longitudinal Framing).ipynb
from math import*
import glob,os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
from matplotlib.font_manager import FontProperties
import matplotlib
import matplotlib.gridspec as gridspec
import time
from gekko import GEKKO
def span_space_stiff(nt,nl,lp,b):
#---------------9. DIMENSIONS OF PANELS AND STIFFENERS------------------------
# Bottom panel
sl = m.Intermediate(1000*b/(nl+1),name='sl')
st = m.Intermediate(1000*lp/(nt+1),name='st')
#if nl < nt transversal arrangements is select by program
#if nt <= nl longitudinal arrangements is select by program
lul = m.Intermediate(st,name='lul')
lut = m.Intermediate(1000*b,name='lut')
return (st, sl, lut, lul)
def stiff_force(lwl,bc,v,beta,mldc,lut,st,lul,sl):
# ---------7. PRESSURE ADJUST FACTORS ACCORDING TO ISO 12215 Part 5-----------
# 7.2.- Design Category Factor (kDC) (Fixed value for Design Category B)
kDC = m.Const(0.8,name='kdcr')
# 7.3.- Dynamic Load Factor (nCG) in "g" units
nCG = (0.32*((lwl/(10*bc))+0.084)*(50-beta)*(((v*bc)**2)/mldc))
if nCG >= 3:
nCG = 0.5*v/(mldc**0.17)
if nCG > 7:
nCG = 7
nCG = m.Const(nCG,name='ncgr')
# 7.5.- Area pressure reduction factor (kAR)
# Design area (AD)
ad_t1 = m.Intermediate(0.33*(lut**2)*1e-6,name='ad_t1')
ad_t2 = m.Intermediate((lut*st)*1e-6,name='ad_t2')
AD_T = m.Intermediate(m.max3(ad_t1,ad_t2),name = 'ad_t')
ad_l1 = m.Intermediate(0.33*(lul**2)*1e-6,name='ad_l1')
ad_l2 = m.Intermediate((lul*sl)*1e-6,name='ad_l2')
AD_L = m.Intermediate(m.max3(ad_l1,ad_l2),name = 'ad_l')
kar_t1 = m.Intermediate((0.1*(mldc**0.15)/(AD_T**0.3)),name='kar_t1')
kar_l1 = m.Intermediate((0.1*(mldc**0.15)/(AD_L**0.3)),name='kar_l1')
# try using m.if3
kar_t2 = m.Intermediate(m.if3(0.25-kar_t1,kar_t1,0.25),name='kar_t2')
kar_l2 = m.Intermediate(m.if3(0.25-kar_l1,kar_l1,0.25),name='kar_l2')
kAR_T = m.Intermediate(m.if3(1-kar_t1,1,kar_t2),name='kar_t')
kAR_L = m.Intermediate(m.if3(1-kar_l1,1,kar_l2),name='kar_l')
#-----------------8. DESIGN PRESSURE ISO 12215 part 5-------------------------
#8.1.3.- Motor craft bottom pressure in planing mode in [kN/m2]
pb_t = m.Intermediate(((0.1*mldc*(1+(kDC**0.5)*nCG))/(lwl*bc))*kAR_T
,name='pb_t')
pb_l = m.Intermediate(((0.1*mldc*(1+(kDC**0.5)*nCG))/(lwl*bc))*kAR_L
,name='pb_l')
#------------11.5 STIMATE DESIGN FORCE AND DESIGN BENDING MOMENT--------------
fd_t = m.Intermediate(5*pb_t*st*lut*1E-4,name='fd_t') #[N]
fd_l = m.Intermediate(5*pb_l*sl*lul*1E-4,name='fd_l') #[N]
md_t = m.Intermediate(83.33*pb_t*st*(lut**2)*1E-9,name='md_t') #[N m]
md_l = m.Intermediate(83.33*pb_l*sl*(lul**2)*1E-9,name='md_l') #[N m]
return pb_t,pb_l,fd_t,fd_l,md_t,md_l
def plate_force(lwl,bc,v,beta,mldc,lut,st,lul,sl,bl,bt):
#---------------------9.DIMENSION OF PANELS AND STIFFENERS--------------------
#9.1.- Dimensions of bottom plating panels
# Large dimension of planel [mm]
l1 = m.max3(sl-bl,st-bt)
# Short dimension of planel [mm]
b1 = m.min3(sl-bl,st-bt)
# ---------7. PRESSURE ADJUST FACTORS ACCORDING TO ISO 12215 Part 5-----------
#7.2.- Design Category Factor (kDC) (Fixed value for Design Category B)
kDC = m.Const(0.8,name='kdcp')
# 7.3.- Dynamic Load Factor (nCG) in "g" units
nCG = 0.32*((lwl/(10*bc))+0.084)*(50-beta)*(((v*bc)**2)/mldc)
if nCG >= 3:
nCG = 0.5*v/(mldc**0.17)
if nCG > 7:
nCG = 7
nCG = m.Const(nCG,name='ncgp')
# 7.5.- Area pressure reduction factor (kAR)
# Design area (AD)
ad1 = m.Intermediate((l1*b1)*1E-6,name='ad1')
ad2 = m.Intermediate(2.5*(b1**2)*1E-6,name='ad2')
AD = m.Intermediate(m.min3(ad1,ad2),name='ad')
# try using m.if3
kar1 = m.Intermediate(0.1*(mldc**0.15)/(AD**0.3),name='kar1')
kar2 = m.Intermediate(m.if3(0.25-kar1,kar1,0.25),name='kar2')
kAR = m.Intermediate(m.if3(1-kar1,1,kar2),name='kar')
#-----------------8. DESIGN PRESSURE ISO 12215 part 5-------------------------
#8.1.3.- Motor craft bottom pressure in planing mode in [kN/m2]
pbp = m.Intermediate(((0.1*mldc*(1+(kDC**0.5)*nCG))/(lwl*bc))*kAR
,name='pbp')
#---------10.1.5 STIMATE FORCE AND BENDING MOMENT ON A PANEL------------------
# Shear Strength aspect ratio factor (kSHC) (Table 12 - ISO 12215 part 5)
kshc1 = m.Intermediate(0.035+0.394*(l1/b1)-0.09*((l1/b1)**2),name='kshc1')
kSHC = m.Intermediate(m.if3(2-(l1/b1),0.5,kshc1),name='kshc')
# Panel Aspec Ratio factor (k2) (Table 5 - ISO 12215 part 5)
a1 = m.Intermediate(0.271*((l1/b1)**2)+0.910*(l1/b1)-0.554,name='a1')
a2 = m.Intermediate(((l1/b1)**2)-0.313*(l1/b1)+1.351,name='a2')
k2 = m.Intermediate(m.if3(2-(l1/b1),0.5,a1/a2),name='k2')
# shear force in middle of "b1" [N/mm]
fdp = m.Intermediate(1*kSHC*pbp*b1*(1E-3),name = 'fdp')
# Bending moment in "b1" direction [N mm/mm]
mdp = m.Intermediate(83.33*2*k2*pbp*(b1**2)*(1E-6),name = 'mdp')
return pbp, fdp, mdp, l1, b1
def bottom_optimize(lp,b,nt,nl,hL,bL,hT,bT,l1,b1,
st,lut,sl,lul,pbt,pbl,fdt,fdl,mdt,mdl,mdp):
#----------STIMATE THE WEIGHT PER AREA, THICKNESS ADN YOUNG'S MODULUS---------
# Bottom Plate
Amat = m.Intermediate(n1*Tmat[0]+n2*Tmat[1]
+n3*Tmat[2]+n4*Tmat[3],name='amat')
Awr = m.Intermediate(n5*Twr[0]+n6*Twr[1],name='awr')
Cmat = m.Intermediate(n1*Wmat[0]+n2*Wmat[1]
+n3*Wmat[2]+n4*Wmat[3],name='cmat')
Cwr = m.Intermediate(n5*Wwr[0]+n6*Wwr[1],name='cwr')
# Young's Modulus of bottom plate laminate [N/mm2]
Ep = m.Intermediate((6400**(4+2*(n1+n2+n3+n4))*(13240**(2*(n5+n6))))
**(1/(4+2*(n1+n2+n3+n4+n5+n6))),name='ep')
# Total laminated thickness [mm]
B_T = m.Intermediate((4*min(Tmat)+2*(Amat+Awr)),name='b_t')
#--------------LONGITUDINAL BOTTOM STIFFENER OPTIMIZATION---------------------
# REQUIERED SHEAR AREA AND SECTIONAL MODULUS
# [cm3] Req. SM in flange of longitudinal stifferners
SM_Req_Long = m.Intermediate((mdl)/(SigUT),name='smrl')
# [cm2] Req. webs shear area of longitudinal stifferners
Awebs_Req_Long = m.Intermediate((fdl)/(100*Tao_MAT),name='awrl')
# [cm2] Req. cores shear area of longitudinal stifferners
Acore_Req_Long = m.Intermediate((fdl)/(100*Tao_Core),name='acrl')
# WEIGHT AND THICKNESS OF BOTTOM LONGITUDINAL
#BT_LS : thickness of bottom longitudinal stiffeners.
BT_LS = m.Intermediate(N_bl1*Tmat[0]+N_bl2*Tmat[1]
+N_bl3*Tmat[2]+N_bl4*Tmat[3],name='bt_ls')
#L_W : weight of bottom longitudinal stiffeners.
L_W = m.Intermediate(N_bl1*Wmat[0]+N_bl2*Wmat[1]
+N_bl3*Wmat[2]+N_bl4*Wmat[3],name='l_w')
# LONGITUDINALS STIFFENERS COMPUTATION
a = m.Intermediate(10+(N_bl1+N_bl2+N_bl3+N_bl4)*15,name='a')
wp = m.Intermediate(20*B_T+bL,name='wp')
EA = m.Intermediate(((Ep*wp*B_T)+(Ec*bL*hL)+
(Er*(a*BT_LS+2*hL*BT_LS+(bL+2*BT_LS)*BT_LS))),
name='EAL')
EAZ = m.Intermediate(((Ep*wp*B_T*(0.5*B_T))+(Ec*bL*hL*(0.5*hL+B_T))+
(Er*((a*BT_LS*(B_T+0.5*BT_LS))+2*hL*BT_LS*(B_T+0.5*hL)
+BT_LS*(bL+2*BT_LS)*(B_T+hL+0.5*BT_LS)))),name='EAZL')
NA = m.Intermediate(EAZ/EA,name='NAL')
Zcri = m.Intermediate(B_T+hL+0.5*BT_LS-NA,name='Zcri')
EAZ2 = m.Intermediate(((Ep*wp*B_T*0.25*(B_T)**2)
+(Ec*((bL*(NA-B_T)*0.25*(NA+B_T)**2)+bL*(hL-NA+B_T)
*(0.25*(NA+hL+B_T)**2)))
+(Er*((2*BT_LS*(NA-B_T)*0.25*(NA+B_T)**2)
+(2*BT_LS*(hL-NA+B_T)*0.25*(NA+hL+B_T)**2)))
+(Er*((a*BT_LS*(B_T+0.5*BT_LS)**2)
+(bL+2*BT_LS)*BT_LS*((B_T+hL+0.5*BT_LS)**2))))
,name='EAZ2L')
Ebh3 = m.Intermediate((1/12)*((Ep*wp*(B_T**3))
+(Ec*bL*(((hL-NA+B_T)**3)+((NA-B_T)**3)))
+(Er*(a*(BT_LS**3)+(bL+2*BT_LS)*(BT_LS**3))))
,name='Ebh3L')
EIna = m.Intermediate(EAZ2+Ebh3-(NA**2)*EA,name='EInaL')
Qc = m.Intermediate(((Ec*0.5*bL*(hL+B_T-NA)**2)
+(Er*(Zcri*BT_LS*(bL+2*BT_LS)+BT_LS*(hL+B_T-NA)**2)))
,name='QcL')
Qw = m.Intermediate(((Ep*(NA-0.5*B_T)*wp*B_T)
+(Er*((NA-B_T-0.5*BT_LS)*a*BT_LS))),name='QwL')
EIminL =m.Intermediate((26*(1**1.5)*pbl*sl*(lul**3)*((1E-7)/0.05)),
name='EIminL')
#DEFINE LONGITUDINAL CONSTRAINTS
#Define longitudinal stiffeners constraints
# [cm3] actual section modulus for longitudinal stiffeners
SML = m.Intermediate((0.001*EIna/(Er*Zcri)),name='sml')
C1L = m.Intermediate((SML/SM_Req_Long),name='c1l')
# [cm2] actual web shear area for longitudinal stiffeners
AWL = m.Intermediate((0.01*EIna*2*BT_LS/Qw),name='awl')
C2L = m.Intermediate((AWL/Awebs_Req_Long),name='c2l')
# [cm2] actual core shear area for longitudinal stiffeners
ACL = m.Intermediate((0.01*EIna*bL/Qc),name='acl')
C3L = m.Intermediate((ACL/Acore_Req_Long),name='c3l')
# [cm2] actual core shear area for longitudinal stiffeners
EIL = m.Intermediate((EIna/EIminL),name='eil')
# realtion for controle the local buckling in webs
C45 = m.Intermediate((hL/BT_LS),name='c45')
# realtion for controle the local buckling in flange
C46 = m.Intermediate((bL/BT_LS),name='c46')
#---------------TRANSVERSAL BOTTOM STIFFENER OPTIMIZATION---------------------
# REQUIERED SHEAR AREA AND SECTIONAL MODULUS
# [cm3] Reg. SM in flange of transversal stiffeners
SM_Req_Trans = m.Intermediate((mdt)/(SigUT),name='smrt')
# [cm2] Req. webs shear area of transversal stiffeners
Awebs_Req_Trans= m.Intermediate((fdt)/(100*Tao_MAT),name='awrt')
# [cm2] Req. cores shear area of transversal stiffeners
Acore_Req_Trans= m.Intermediate((fdt)/(100*Tao_Core),name='acrt')
# WEIGHT AND THICKNESS OF BOTTOM TRANSVERSAL
#BT_TS : thickness of bottom transversal stiffeners.
BT_TS = m.Intermediate(N_bt1*Tmat[0]+N_bt2*Tmat[1]
+N_bt3*Tmat[2]+N_bt4*Tmat[3],name='bt_ts')
#T_W : weight of bottom transversal stiffeners.
T_W = m.Intermediate(N_bt1*Wmat[0]+N_bt2*Wmat[1]
+N_bt3*Wmat[2]+N_bt4*Wmat[3],name='t_w')
# TRANSVERSAL STIFFENERS COMPUTATION
aT =m.Intermediate((10+(N_bt1+N_bt2+N_bt3+N_bt4)*15),name='aT')
wpT =m.Intermediate(20*B_T+bT,name='wpT')
EAT =m.Intermediate((Ep*(wpT*B_T)+(Ec*bT*hT)+
(Er*(aT*BT_TS+2*hT*BT_TS+(bT+2*BT_TS)*BT_TS))),
name='EAT')
EAZT =m.Intermediate(((Ep*wpT*B_T*(0.5*B_T))+(Ec*bT*hT*(0.5*hT+B_T))+
(Er*((aT*BT_TS*(B_T+0.5*BT_TS))+
2*hT*BT_TS*(B_T+0.5*hT)+
BT_TS*(bT+2*BT_TS)*(B_T+hT+0.5*BT_TS)))),
name='EAZT')
NAT = m.Intermediate(EAZT/EAT,name='NAT')
ZcriT =m.Intermediate(B_T+hT+0.5*BT_TS-NAT,name='ZcriT')
EAZ2T =m.Intermediate(((Ep*wpT*B_T*0.25*(B_T**2))+
(Ec*(bT*(NAT-B_T)*0.25*((NAT+B_T)**2)+
(bT*(hT-NAT+B_T)*(0.25*(NAT+hT+B_T)**2))))+
(Er*((2*BT_TS*(NAT-B_T)*0.25*(NAT+B_T)**2)+
(2*BT_TS*(hT-NAT+B_T)*0.25*(NAT+hT+B_T)**2)))+
(Er*((aT*BT_TS*(B_T+0.5*BT_TS)**2)+
(bT+2*BT_TS)*BT_TS*((B_T+hT+0.5*BT_TS)**2)))),
name='EAZ2T')
Ebh3T =m.Intermediate((1/12)*((Ep*wpT*(B_T**3))+
(Ec*bT*(((hT-NAT+B_T)**3)+((NAT-B_T)**3)))+
(Er*2*BT_TS*(((hT-NAT+B_T)**3)+
((NAT-B_T)**3)))+
(Er*(aT*(BT_TS**3)+(bT+2*BT_TS)*(BT_TS**3)))),
name='Ebh3T')
EInaT =m.Intermediate(EAZ2T+Ebh3T-(NAT**2)*EAT,name='EInaT')
Qc_T =m.Intermediate(((Ec*0.5*bT*(hT+B_T-NAT)**2)
+(Er*(ZcriT*BT_TS*(bT+2*BT_TS)+
BT_TS*(hT+B_T-NAT)**2))),name='QcT')
Qw_T =m.Intermediate(((Ep*(NAT-0.5*B_T)*wpT*B_T)+
(Er*((NAT-B_T-0.5*BT_TS)*aT*BT_TS))),name='QwT')
EIminT =m.Intermediate((26*(1**1.5)*pbt*st*(lut**3)*((1E-7)/0.05)),
name='EIminT')
# DEFINE TRANSVERSAL CONSTRAINTS
# [cm3] actual section modulus for transversal stiffeners
SMT = m.Intermediate(0.001*EInaT/(Er*ZcriT),name='smt')
C1T = m.Intermediate((SMT/SM_Req_Trans),name='c1t')
# [cm2] actual web shear area for transversal stiffeners
AWT = m.Intermediate((0.01*EInaT*2*BT_TS/Qw_T),name='awt')
C2T = m.Intermediate((AWT/Awebs_Req_Trans),name='c2t')
# [cm2] actual core shear area for transversal stiffeners
ACT = m.Intermediate((0.01*EInaT*bT/Qc_T),name='act')
C3T = m.Intermediate((ACT/Acore_Req_Trans),name='c3t')
# relation of stiffeness for transversal stiffeners
EIT = m.Intermediate((EInaT/EIminT),name='eit')
# relation for controle the local buckling in webs
C47 = m.Intermediate((hT/BT_TS),name='c47')
# relation for controle the local buckling in flange
C48 = m.Intermediate((bT/BT_TS),name='c48')
#----------------------TOTAL STIFFENERS WEIGHT--------------------------------
L_Weight = m.Intermediate(((L_W*nl*lul/300)*((2*hL+bL+a)/1000)+
(hL*bL*lul*nl*dcore/(1000**3))),name='LS_weight')
T_Weight = m.Intermediate(((T_W*nt*lut/0.3)*((2*hT+bT+aT)/1000)
+(hT*bT*nt*lut*dcore/1000**2))
,name='TS_weight')
#weight per length of mat in longitudinal stiffeners.
lwl = m.Intermediate((L_W/0.3)*((2*hL+bL+a)/1000),name='lwlayer')
#weight per length of core in longitudinal stiffeners.
lwc = m.Intermediate(hL*bL*dcore/1000**2,name='lwcore')
#weight per length of mat in transversal stiffeners
twl = m.Intermediate((T_W/0.3)*((2*hT+bT+aT)/1000),name='twlayer')
#weight per length of core in transversal stiffeners
twc = m.Intermediate(hT*bT*dcore/1000**2,name='twcore')
#---------------------------PLATE OPTIMIZATIONS-------------------------------
# PLATE COMPUTATION
# Sectional Inercia
I = m.Intermediate((((4*min(Tmat)+2*(Amat+Awr))**3)/12),name='I')
# Distance to critical fiber
Zcri = m.Intermediate((Amat+Awr)+1.5*min(Tmat),name='Zcrip')
# Acting tension stress
Sact = m.Intermediate((mdp/(Ep*I/(Zcri*6400))),name='Sact')
# Total Panel dimension in squared meters
Panel_Dim = lp*b
# Local panel dimension in squared meters
LPanel_Dim = m.Intermediate(l1*b1/(1000**2),name='Lpanel')
# CONSTRAINTS
#Sub-Laminate
C1P = m.Intermediate(n1+n2+n3+n4-n5-n6,name='c1p')
# Min. Weight per area [kg/m2]
minWeight = m.Const((0.0675+0.0135*v*0.0675*(mldc**0.33)),name='minW')
bweight = m.Intermediate((4*min(Wmat)/0.3)+
2*((Cmat/0.3)+(Cwr/0.48)),name='bweight')
C2P = m.Intermediate((minWeight/bweight),name='c2p')
#Norma stress acting on critical fiber
C3P = m.Intermediate((SigUT/Sact),name='c3p')
#---------OBJECTIVE FUNCTION FOR BOTTOM PANEL (panel weight in Kg )-----------
PANEL_WEIGHT = m.Intermediate(Panel_Dim*bweight,name='pw')
STIFF_WEIGHT = m.Intermediate((L_Weight+T_Weight),name='sw')
BOTTOM_WEIGHTS = PANEL_WEIGHT+STIFF_WEIGHT
return (BOTTOM_WEIGHTS,PANEL_WEIGHT,STIFF_WEIGHT,
C1L,C2L,C3L,EIL,C45,C46,C1T,C2T,C3T,EIT,C47,C48,C1P,C2P,C3P,
B_T,Ep,twl,twc,lwl,lwc)
#=========================DEFINE SHIP PARAMETERS================================
LH = 12.0 # [m] Hull length
Lwl = 10.74 # [m] Legnth on design water line
BC = 2.93 # [m] Chine beam in meters
beta = 13 # [°] Deadrise angle
V = 25 # [kn] Maximun velocity
mLDC = 8864 # [kg] Design Displacement
T = 0.51 # [m] Draft
D = 1.5 # [m] Depth
#=============================FUNCTION PARAMETERS===============================
#DEFINE FUNCTION PARAMETERS
Lp = 4 # 4x1.5 [m] global bottom panel dimension in "x" direction (Aft/Fwr)
B = 1.5 # 4x1.5 [m] global bottom panel dimension in "y" direction (Br/Stbr)
Bs = 1.2 # 4x1.2 [m] global side panel dimension in "y" direction
# it is necessary change the name of variable for functions COMPUTATION
lh,lwl,bc,v,mldc,lp,b,bside,d,t = LH,Lwl,BC,V,mLDC,Lp,B,Bs,D,T
#Permisibles stress, (Reverse Safe Factor = 0.4) (ISO 12215 part 5 Annex C)
SigUT = 0.4*85 # N/mm2 ultimate tension stress MAT
Taoi_MAT = 0.4*17.25 # N/mm2 ultimate intralaminar stress MAT
Taoi_WR = 0.4*14.1 # N/mm2 ultimate intralaminar stress WR
Tao_MAT = 0.4*62 # N/mm2 ultimate shear strength stress MAT
Tao_Core = 0.4*16.7 # N/mm2 ultimate shear strength stress CORE
# Young's Modulus
Ec = 17240 # N/mm2 young's modulus for Core material (CHANUL)
Er = 6400 # N/mm2 young's modulus for MAT fibers
#Material properties (ISO 12215 part 5 Annex C)
dres = 1.2 # g/cm3 Polyester resine density
dmat = 2.56 # g/cm3 Mat ply density
dwr = 2.56 # g/cm3 Waven Roven ply density
Wmat =[0.3,0.35,0.45,0.6] # weigth per area (kg/m2)
# of Mat plies (300,350,450,600) grams.
Wwr =[0.6,0.8] # weigth per area (kg/m2) of WR plies (600,800) grams
dkeel = 870 # Kg/m3 Keel Core density (chanul)
dcore = 870 # Kg/m3 Stiffeners Core density (chanul)
#Thickness of each layer
Tmat = np.zeros(len(Wmat))
Twr = np.zeros(len(Wwr))
for i in range(len(Wmat)):
Tmat[i] = (Wmat[i]/(dres*dmat))*((dmat/0.3)-(dmat-dres))
for i in range(len(Wwr)):
Twr[i] = (Wwr[i]/(dres*dwr))*((dwr/0.48)-(dwr-dres))
#======================OPTIMIZATION OF BOTTOM PANELS============================
#Initial point for bottom optimization
n1, n2, n3, n4, n5, n6 = 10,10,2,2,10,1
NL, NT =1,1
N_bl1, N_bl2, N_bl3, N_bl4 = 0,5,1,2
N_bt1, N_bt2, N_bt3, N_bt4 = 2,0,1,2
hL ,bL ,hT ,bT = 150,100,60,40
#Vectors with initial points
Nmat = [n1,n2,n3,n4] # number of MAT layers
Nwr = [n5,n6] # number of WR layers
BL_Mat = [N_bl1, N_bl2, N_bl3, N_bl4] # longitudinal stiffeners laminated
BT_Mat = [N_bt1, N_bt2, N_bt3, N_bt4] # transversal stiffeners laminated
DIM = [hL ,bL ,hT ,bT] # stiffeners dimensions
na= [n1,n2,n3,n4,n5,n6] # number of layer in bottom panel
NR = [NT,NL] # number of bottom stiffeners
#====================== GEKKO FIRST FASE ITERACTIONS ===========================
print(74*"=")
print(74*"=")
print(28*"=","START OPTIMIZATION",28*"=")
start_time = time.time()
print(74*"=")
from gekko import GEKKO
m = GEKKO(remote=False) # Initialize gekko
m.options.IMODE = 3 # solution mode 2 or 3
m.options.SOLVER= 1 # (1=APOPT, 3=IPOPT) APOPT is an MINLP solver
# optional solver settings with APOPT
m.solver_options = ['minlp_maximum_iterations 6000', \
# minlp iterations with integer solution
'minlp_max_iter_with_int_sol 6000', \
# treat minlp as nlp
'minlp_as_nlp 0', \
# nlp sub-problem max iterations
'nlp_maximum_iterations 5000', \
# 1 = depth first, 2 = breadth first , 3 , 4
'minlp_branch_method 1', \
#maximum value to be considered as an integer
'minlp_integer_max 2.0e9'\
# maximum deviation from whole number
'minlp_integer_tol 1e-5', \
# covergence tolerance
'minlp_gap_tol 0.001'\
#objective_convergence_tolerance 1.0e-6
'constraint_convergence_tolerance 1.0e-6'\
#convergence tolerance for the constraints
'objective_convergence_tolerance 1.0e-6']
#======================= PRINT INITIAL POINTS. =================================
print(8*"..","Initial point for bottom optimization",9*"..")
print("Number of Stiffeners [NT,NL] : ",NR)
print("")
print("Bottom Layers [n1,n2,n3,n4,n5,n6] : ",na)
print("")
print("Long. Stiffeners layers [bln1,bln2,bln3,bln4] : ",BL_Mat)
print("")
print("Trans. Stiffeners layer [btn1,btn2,btn3,btn4] : ",BT_Mat)
print("")
print("Stiffeners dimension [hL,bL,hT,bT][mm] : ",DIM)
print(37*"..")
#============================ DESIGN VARIABLES.=================================
low=1
n1 = m.Var(value=na[0],lb=low,ub=100,integer=True,name='n1')
n2 = m.Var(value=na[1],lb=low,ub=100,integer=True,name='n2')
n3 = m.Var(value=na[2],lb=low,ub=100,integer=True,name='n3')
n4 = m.Var(value=na[3],lb=low,ub=100,integer=True,name='n4')
n5 = m.Var(value=na[4],lb=low,ub=100,integer=True,name='n5')
n6 = m.Var(value=na[5],lb=low,ub=100,integer=True,name='n6')
NT = m.Var(value=NR[0],lb=low,ub=10,integer=True,name='NT')
NL = m.Var(value=NR[0],lb=low,ub=3 ,integer=True,name='NL')
N_bl1 = m.Var(value=BL_Mat[0],lb=low,ub=100,integer=True,name='N_bl1')
N_bl2 = m.Var(value=BL_Mat[1],lb=low,ub=100,integer=True,name='N_bl2')
N_bl3 = m.Var(value=BL_Mat[2],lb=low,ub=100,integer=True,name='N_bl3')
N_bl4 = m.Var(value=BL_Mat[3],lb=low,ub=100,integer=True,name='N_bl4')
N_bt1 = m.Var(value=BT_Mat[0],lb=low,ub=100,integer=True,name='N_bt1')
N_bt2 = m.Var(value=BT_Mat[1],lb=low,ub=100,integer=True,name='N_bt2')
N_bt3 = m.Var(value=BT_Mat[2],lb=low,ub=100,integer=True,name='N_bt3')
N_bt4 = m.Var(value=BT_Mat[3],lb=low,ub=100,integer=True,name='N_bt4')
hL = m.Var(value=DIM[0],lb=40,ub=800,integer=False,name='hl')
bL = m.Var(value=DIM[1],lb=40,ub=800,integer=False,name='bl')
hT = m.Var(value=DIM[2],lb=40,ub=800,integer=False,name='ht')
bT = m.Var(value=DIM[3],lb=40,ub=800,integer=False,name='bt')
#============================GEKKO COMPUTATION.=================================
#Stimation of space and span stiffeners
sT, sL, luT, luL = span_space_stiff(NT,NL,Lp,B)
#Stimation of stiffeners forces
pb_T,pb_L,fd_T,fd_L,md_T,md_L = stiff_force(Lwl,BC,V,beta,mLDC,luT,sT,luL,sL)
#Stimation of plate bending moment and shear force
pb_p,fd_p,md_p,l_1,b_1 = plate_force(Lwl,BC,V,beta,mLDC,luT,sT,luL,sL,bL,bT)
#Plate optimiation.
(BOTTOM_WEIGHTS,PANEL_WEIGHT,STIFF_WEIGHT,
C1L,C2L,C3L,EIL,C45,C46,C1T,C2T,C3T,EIT, C47,C48,C1P,C2P,C3P,
B_T,EP,tw_l,tw_c,lw_l,lw_c)=bottom_optimize(Lp,B,NT,NL,hL,bL,hT,bT,l_1,b_1,sT,
luT,sL,luL,pb_T,pb_L,fd_T,fd_L,
md_T,md_L,md_p)
#=============================CONSTRAINTS.======================================
# Panel constraints.
# Same number of Mat layers and Wr layers into sub laminate
CON1 = m.Equation(C1P == 0)
# Minimun weight for impact loads
CON2 = m.Equation(C2P -1 >= 0)
#Define Objective function for panel weight minimization
Objective = m.Obj(BOTTOM_WEIGHTS)
m.open_folder()
output = m.solve(disp=True) # Solver
Try setting a higher diagnostic level such as
DIAGLEVEL=4:In the run folder, you'll see a couple new files that are helpful in diagnosing the problems with the equations. The first is
model_init_values.txtthat shows the initial guesses for the values with the intermediates calculated based on the guess values that you provided. You'll see a number ofNaNorInfvalues that likely mean divide by zero.Another helpful file is
model_init_res.txtthat shows the equation residuals (f(x)==g(x)residual isf(x)-g(x)). These are some of the equation residuals:I recommend that you examine each of the equations and provide an initial guess that produces a finite values. Also, it is a good idea to protect against divide by zero. For example, you can switch an Intermediate equation like
y = m.Intermediate(1/x)to a variable foryand reformulate the equation asm.Equation(y*x==1). I'll be glad to provide additional help on this problem - it looks like a very nice application.