【发布时间】：2020-12-28 11:54:02
【问题描述】：
我正在尝试使用 MINLP 和 GEKKO 包开发玻璃增强塑料 (GRP) 高速工艺的结构重量优化。我需要将许多设计变量视为整数，而将其他设计变量视为实数和连续变量。为了评估约束，我遵循 ISO 12215 第 5 部分的要求。该规则确定了工艺的最低结构要求。问题是我需要使用许多条件，例如（if、max 或 min），我读到 GEKKO 改变了这些类型的条件语句，因为传统的条件不是连续的。因此，我使用m.if3、m.max3 或m.min3。但是，当我运行程序时，这些 GEKKO 条件语句不起作用。许多运算（中间体）取无穷大值，因为程序除以零并且永远不会收敛到解。有什么建议可以理解求解器并获得成功的解决方案？
我的问题相对较短，因为有 20 个设计变量和大约 15 个约束，没有时间依赖性。
# 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
【问题讨论】：
-
感谢您的提问，米格尔。我将您的问题编辑为更具体。这里有一些建议可以帮助你：（1）查看run文件夹中的infeasibilties.txt文件。您可以包含
m.Var(name='myVar')以使该文件更易于阅读 (2) 通过将x/y==1重新公式化为x==y来避免除以零的方程式（您可能需要先使用m.Equation而不是m.Intermediate(3)通过将设计变量从m.Var切换到m.Param来找到可行的初始猜测值。此问题已被版主关闭。如果您仍需要帮助，请再次发布。 -
这是一个相关问题，将
y==1/x重新表述为y*x==1有助于解决问题：stackoverflow.com/questions/64491175/… 有时Variable+Equation比Intermediate更容易解决因为它避免被零除。另外，检查if3()或max3()函数是否可以被变量约束替换。这是一个相关问题，它没有您的复杂，但具有一些相同的设计特征。 apmonitor.com/me575/index.php/Main/TubularColumn
标签： python-3.x optimization gekko