Miguel
Reputation: 21
Python Mixed Integer Nonlinear Optimization using GEKKO package
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
Upvotes: 1
Views: 800
Answers (1)
John Hedengren
Reputation: 14346
Try setting a higher diagnostic level such as DIAGLEVEL=4:
m.open_folder()
m.options.DIAGLEVEL=4
output = m.solve(disp=True) # Solver
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.txt that shows the initial guesses for the values with the intermediates calculated based on the guess values that you provided. You'll see a number of NaN or Inf values that likely mean divide by zero.
172 +Inf kar1
173 0.00000000 kar2
174 0.00000000 kar
175 0.00000000 pbp
176 NaN kshc1
177 0.00000000 kshc
178 NaN a1
179 NaN a2
180 0.00000000 k2
181 0.00000000 fdp
182 0.00000000 mdp
.
.
.
195 1311.11111111 wp
196 917050883.35535383 eal
197 57906932528.31179047 eazl
198 63.14473229 nal
199 151.38044132 zcri
200 ******************** eaz2l
201 ******************** ebh3l
202 ******************** einal
203 20726655556.62884521 qcl
204 20720536220.47989273 qwl
205 0.00000000 eiminl
206 3025.53188955 sml
207 +Inf c1l
208 22.46255737 awl
209 +Inf c2l
210 141.42371638 acl
211 +Inf c3l
212 +Inf eil
213 18.89350536 c45
Another helpful file is model_init_res.txt that shows the equation residuals (f(x)==g(x) residual is f(x)-g(x)). These are some of the equation residuals:
6 116 -1.32000000 v24-(((((1-int_v23))*(ad_l1))+((int_v23)*(ad_l2))))
7 117 +Inf 0-((((1-int_v25))*((0.25-kar_t1))))-slk_7
8 118 -Inf ((int_v25)*((0.25-kar_t1)))-(0)-slk_8
13 123 +Inf 0-((((1-int_v29))*((1-kar_t1))))-slk_13
14 124 -Inf ((int_v29)*((1-kar_t1)))-(0)-slk_14
15 125 -0.99000000 v30-(((((1-int_v29))*(1))+((int_v29)*(kar_t2))))
16 126 +Inf 0-((((1-int_v31))*((1-kar_l1))))-slk_16
17 127 -Inf ((int_v31)*((1-kar_l1)))-(0)-slk_17
18 128 -0.99000000 v32-(((((1-int_v31))*(1))+((int_v31)*(kar_l2))))
19 129 -1310.00000000 0-((((1-int_v33))*(((st-bt)-(sl-bl)))))-slk_19
26 136 0.00000000 0-(((int_v37)*((ad2-ad1))))-slk_26
27 137 0.00000000 v38-(((((1-int_v37))*(ad1))+((int_v37)*(ad2))))
28 138 +Inf 0-((((1-int_v39))*((0.25-kar1))))-slk_28
29 139 -Inf ((int_v39)*((0.25-kar1)))-(0)-slk_29
30 140 -Inf v40-(((((1-int_v39))*(kar1))+((int_v39)*(0.25))))
33 143 -0.99000000 v42-(((((1-int_v41))*(1))+((int_v41)*(kar2))))
34 144 NaN 0-((((1-int_v43))*((2-((v34)/(v36))))))-slk_34
35 145 NaN ((int_v43)*((2-((v34)/(v36)))))-(0)-slk_35
40 150 13.00000000 c1p-(0)
41 151 -0.99414611 (c2p-1)-(0)-slk_41
42 152 6464.36333333 (pw+sw)
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 for y and reformulate the equation as m.Equation(y*x==1). I'll be glad to provide additional help on this problem - it looks like a very nice application.
Upvotes: 0
Related Questions
- Mixed integer nonlinear programming with gekko python
- mixed integer programming optimization
- Formulate a nonlinear optimization problem in Python (gekko) and with m.Equations
- Gekko optimisation error with a nonlinear function and 2 inputs
- Gekko : MINLP options not accepted by gekko
- Issue with integer optimization with Gekko
- optimization using package
- Calling Gekko solve gives TypeError: object of type 'int' has no len()
- Gekko optimization
- Gekko non linear mixed integer programming. Problem setting up objective function