📄 form1.frm
字号:
VERSION 5.00
Begin VB.Form Form1
Caption = "Form1"
ClientHeight = 3090
ClientLeft = 60
ClientTop = 450
ClientWidth = 4680
LinkTopic = "Form1"
ScaleHeight = 3090
ScaleWidth = 4680
StartUpPosition = 3 '窗口缺省
Begin VB.CommandButton Command1
Caption = "Command1"
Height = 615
Left = 2520
TabIndex = 0
Top = 1560
Width = 1815
End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Option Base 1
Dim x1a() As Variant, x2a() As Variant, x3a() As Variant
Dim x1b() As Variant, x2b() As Variant, x3b() As Variant
Dim xa(DataNum, 3) As Double, xb(DataNum, 3) As Double
Private Sub Command1_Click()
Dim a() As Double, b() As Double, VaryNum%, Popsize%, MaxGeneration%, ObjectValue As Integer, M As Double
'a() 变量的取值下限 b() 变量的取值上限
'VaryNum 变量数 MaxGeneration 最大迭代数目
'Popsize 微粒群规模
VaryNum = 6 '微粒群算法中变量总数
ReDim a(VaryNum) As Double, b(VaryNum) As Double
For i = 1 To VaryNum
a(i) = -4
b(i) = 4
Next i
Popsize = 150 '微粒群算法中个体总数
MaxGeneration = 5000 '微粒群算法运行代数
ObjectValue = 0 '判断目标函数取最小值还是最大值;当取最小值时取0,取最大值时取1
M = Pso(a(), b(), VaryNum, Popsize, MaxGeneration, ObjectValue)
End Sub
Private Sub Form_Load() '数据录入
'##########数据来源 /徐桦等,水－正已烷－乙醇体系的液液平衡研究.,20组308.15K时LLE数据/############
x1a = Array(0.9817, 0.9781, 0.9762, 0.9677, 0.9559, 0.9432, 0.9259, 0.8827, 0.8549, 0.8239, 0.8009, 0.7759, 0.7524, 0.629, 0.5013, 0.3429, 0.2228, 0.1574, 0.125, 0.0784)
'水相水含量x1a
x2a = Array(0.018, 0.0212, 0.0225, 0.031, 0.0425, 0.055, 0.0723, 0.1151, 0.1398, 0.171, 0.1929, 0.2102, 0.2306, 0.35, 0.467, 0.5999, 0.6825, 0.7054, 0.58, 0.49)
'水相乙醇含量x2a
x3a = Array(0.0003, 0.0007, 0.0013, 0.0013, 0.0016, 0.0018, 0.0018, 0.0022, 0.0053, 0.0051, 0.0062, 0.0139, 0.017, 0.021, 0.0317, 0.0572, 0.0947, 0.1372, 0.295, 0.4316)
'水相正己烷含量x3a
x1b = Array(0.0006, 0.0002, 0.0003, 0.0007, 0.0004, 0.0005, 0.0006, 0.0019, 0.0029, 0.0044, 0.0052, 0.0014, 0.0007, 0.0045, 0.006, 0.0064, 0.0062, 0.0028, 0.0057, 0.0064)
'油相水含量x1b
x2b = Array(0.0008, 0.0014, 0.003, 0.0023, 0.0099, 0.0032, 0.0147, 0.0155, 0.0165, 0.0182, 0.0222, 0.0275, 0.0314, 0.033, 0.0342, 0.0358, 0.0364, 0.0372, 0.0345, 0.0392)
'油相乙醇含量x2b
x3b = Array(0.9986, 0.9984, 0.9967, 0.997, 0.9897, 0.9963, 0.9847, 0.9826, 0.9806, 0.9774, 0.9726, 0.9711, 0.9679, 0.9625, 0.9598, 0.9578, 0.9574, 0.96, 0.9598, 0.9544)
'油相正己烷含量x3b
'##########数据来源 /醋酸甲酯（1）－甲醇（2）－水（3）三元液液相平衡数据,293.15K时LLE数据/############
'x1a = Array(0, 0.0771, 0.1067, 0.1315, 0.1584, 0.2063)
'x2a = Array(0, 0.002, 0.0359, 0.0533, 0.0628, 0.0706)
'x3a = Array(0, 0.9209, 0.8574, 0.8152, 0.7788, 0.7231)
'x1b = Array(0, 0.7515, 0.6851, 0.6335, 0.4966, 0.3744)
'x2b = Array(0, 0.0062, 0.035, 0.0441, 0.0717, 0.0724)
'x3b = Array(0, 0.2433, 0.2799, 0.3224, 0.4317, 0.5532)
For i = 1 To DataNum
xa(i, 1) = x1a(i)
xa(i, 2) = x2a(i)
xa(i, 3) = x3a(i)
xb(i, 1) = x1b(i)
xb(i, 2) = x2b(i)
xb(i, 3) = x3b(i)
Next i
End Sub
Public Function FunctionObject(Vary() As Double, Generation As Integer) As Double
Dim Dimension As Integer, a() As Double, i%, j%
Dimension = UBound(Vary)
ReDim a(Dimension) As Double
For i = 1 To Dimension
a(i) = Vary(i)
Next i
'以下目标函数变时需要修改
Dim fValue(DataNum) As Double
Dim G(6) As Double, M#
Dim ra(DataNum, 3) As Double, rb(DataNum, 3) As Double '定义活度系数
Dim Y11(DataNum) As Double, Y12(DataNum) As Double, Y13(DataNum) As Double, Y14(DataNum) As Double, Y15(DataNum) As Double, Y16(DataNum) As Double
Dim Y21(DataNum) As Double, Y22(DataNum) As Double, Y23(DataNum) As Double, Y24(DataNum) As Double, Y25(DataNum) As Double, Y26(DataNum) As Double
Dim Y31(DataNum) As Double, Y32(DataNum) As Double, Y33(DataNum) As Double, Y34(DataNum) As Double, Y35(DataNum) As Double, Y36(DataNum) As Double
Dim KReal(DataNum, 3) As Double, KCal(DataNum, 3) As Double
'Nrtl方程参数求解
M = 0.2 'm 代表NRTL方程第三参数，根据文献取值为0.3
G(1) = Exp(-M * a(1)) 'G(i) 代表NRTL方程参数
G(2) = Exp(-M * a(2)) 'a(i) 代表NRTL方程参数τ
G(3) = Exp(-M * a(3)) 'a(1)代表τ12,a(2)代表τ13
G(4) = Exp(-M * a(4)) 'a(3)代表τ21,a(4)代表τ23
G(5) = Exp(-M * a(5)) 'a(5)代表τ31,a(6)代表τ32
G(6) = Exp(-M * a(6))
' x1a代表α相(水相)中水含量
' x2a代表α相(水相)中乙醇含量
' x3a代表α相(水相)中正己烷含量
'########求lnr1(求α、β相中活度系数lnr1)###########
For i = 1 To DataNum
Y11(i) = (x2b(i) * a(3) * G(3) + x3b(i) * a(5) * G(5)) * (x2b(i) * G(3) + x3b(i) * G(5)) / (x1b(i) + x2b(i) * G(3) + x3b(i) * G(5)) ^ 2
Y12(i) = G(1) * x2b(i) * (a(1) * x2b(i) + a(1) * x3b(i) * G(6) - x3b(i) * a(6) * G(6)) / (x1b(i) * G(1) + x2b(i) + x3b(i) * G(6)) ^ 2
Y13(i) = G(2) * x3b(i) * (x2b(i) * a(2) * G(2) + x3b(i) * a(2) - x2b(i) * a(4) * G(4)) / (x1b(i) * G(2) + x2b(i) * G(4) + x3b(i)) ^ 2
Y14(i) = (x2a(i) * a(3) * G(3) + x3a(i) * a(5) * G(5)) * (x2a(i) * G(3) + x3a(i) * G(5)) / (x1a(i) + x2a(i) * G(3) + x3a(i) * G(5)) ^ 2
Y15(i) = G(1) * x2a(i) * (a(1) * x2a(i) + a(1) * x3a(i) * G(6) - x3a(i) * a(6) * G(6)) / (x1a(i) * G(1) + x2a(i) + x3a(i) * G(6)) ^ 2
Y16(i) = G(2) * x3a(i) * (x2a(i) * a(2) * G(2) + x3a(i) * a(2) - x2a(i) * a(4) * G(4)) / (x1a(i) * G(2) + x2a(i) * G(4) + x3a(i)) ^ 2
Next i
For i = 1 To DataNum
rb(i, 1) = Exp(Y11(i) + Y12(i) + Y13(i))
ra(i, 1) = Exp(Y14(i) + Y15(i) + Y16(i))
Next i
'########求lnr1(求α、β相中活度系数lnr1)#############
'########求lnr2(求α、β相中活度系数lnr2)##########
For i = 1 To DataNum
Y21(i) = (x1b(i) * G(1) * a(1) + x3b(i) * a(6) * G(6)) * (x1b(i) * G(1) + x3b(i) * G(6)) / (x1b(i) * G(1) + x2b(i) + x3b(i) * G(6)) ^ 2
Y22(i) = x1b(i) * G(1) * (x1b(i) * a(3) + x3b(i) * a(3) * G(5) - x3b(i) * a(5) * G(5)) / (x1b(i) + x2b(i) * G(3) + x3b(i) * G(5)) ^ 2
Y23(i) = x3b(i) * G(4) * (a(4) * x3b(i) + a(4) * x1b(i) * G(2) - x1b(i) * G(2) * a(2)) / (x1b(i) * G(2) + x2b(i) * G(4) + x3b(i)) ^ 2
Y24(i) = (x1a(i) * a(1) * G(1) + x3a(i) * a(6) * G(6)) * (x1a(i) * G(1) + x3a(i) * G(6)) / (x1a(i) * G(1) + x2a(i) + x3a(i) + G(6)) ^ 2
Y25(i) = x1a(i) * G(1) * (a(3) * x1a(i) + x3a(i) * a(3) * G(5) - x3a(i) * a(5) * G(5)) / (x1a(i) + x2a(i) * G(3) + x3a(i) * G(5)) ^ 2
Y26(i) = x3a(i) * G(4) * (a(4) * x3a(i) + a(4) * x1a(i) * G(2) - x1a(i) * a(2) * G(2)) / (x1a(i) * G(2) + x2a(i) * G(4) + x3a(i)) ^ 2
'表达式2 lnγ2I=Y21+Y22+Y23 lnγ2II=Y24+Y25+Y26 Y2=ln(γ2I/γ2II)
Next i
For i = 1 To DataNum
rb(i, 2) = Exp(Y21(i) + Y22(i) + Y23(i))
ra(i, 2) = Exp(Y24(i) + Y25(i) + Y26(i))
Next i
'########求lnr2(求α、β相中活度系数lnr2)##########
'########求lnr3(求α、β相中活度系数lnr3)##########
For i = 1 To DataNum
Y31(i) = (x1b(i) * a(2) * G(2) + x2b(i) * a(4) * G(4)) * (x1b(i) * G(2) + x2b(i) * G(4)) / (x1b(i) * G(2) + x2b(i) * G(4) + x3b(i)) ^ 2
Y32(i) = x1b(i) * G(5) * (a(5) * x1b(i) + a(5) * x2b(i) * G(3) - x2b(i) * G(3) * a(3)) / (x1b(i) + x2b(i) * G(3) + x3b(i) * G(5)) ^ 2
Y33(i) = x2b(i) * G(6) * (x2b(i) * a(6) + x1b(i) * a(6) * G(6) - x1b(i) * G(1) * a(1)) / (x1b(i) * G(1) + x2b(i) + x3b(i) * G(6)) ^ 2
Y34(i) = (x1a(i) * a(2) * G(2) + x2a(i) * G(4) * a(4)) * (x1a(i) * G(2) + x2a(i) * G(4)) / (x1a(i) * G(2) + x2a(i) * G(4) + x3a(i)) ^ 2
Y35(i) = x1a(i) * G(5) * (a(5) * x1a(i) + a(5) * x2a(i) * G(3) - x2a(i) * a(3) * G(3)) / (x1a(i) + x2a(i) * G(3) + x3a(i) * G(5)) ^ 2
Y36(i) = x2a(i) * G(6) * (x2a(i) * a(6) + x1a(i) * a(6) * G(6) - x1a(i) * a(1) * G(1)) / (x1a(i) * G(1) + x2a(i) + x3a(i) * G(6)) ^ 2
'表达式3 lnγ3I=Y31+Y32+Y33 lnγ3II=Y34+Y35+Y36 Y3=ln(γ3I/γ3II)
Next i
For i = 1 To DataNum
rb(i, 3) = Exp(Y31(i) + Y32(i) + Y33(i))
ra(i, 3) = Exp(Y34(i) + Y35(i) + Y36(i))
Next i
'########求lnr3(求α、β相中活度系数lnr2)##########
For i = 1 To DataNum
For j = 1 To 3
KReal(i, j) = xa(i, j) / xb(i, j)
KCal(i, j) = rb(i, j) / ra(i, j)
Next j
Next i
FunctionObject = 0 '另初值为零，得注意这里
For i = 1 To DataNum
For j = 1 To 3
FunctionObject = FunctionObject + (Abs((KReal(i, j) - KCal(i, j)) / KReal(i, j))) ^ 2 '目标函数 值越小越好
Next j
Next i
FunctionObject = FunctionObject / (3 * DataNum)
'以上目标函数变时需要修改
End Function
Public Function Pso(amin() As Double, amax() As Double, VaryNum As Integer, Popsize As Integer, MaxGeneration As Integer, ObjectValue As Integer)
'amin 变量的取值下限 amax 变量的取值上限
'VaryNum 变量数 MaxGeneration 最大迭代数目
'Popsize 微粒群规模
Dim RealValue(DataNum) As Double
Dim i%, j%, Generation%, PNum%
Dim c1#, c2#, Rnd1#, Rnd2#, w#
'##########重定义数组变量##############
Dim InBird() As Individual, IndividualBest() As Individual, GlobalBest As Individual
ReDim GlobalBest.v(VaryNum) As Double
ReDim GlobalBest.x(VaryNum) As Double
ReDim InBird(MaxGeneration, Popsize) As Individual, IndividualBest(Popsize) As Individual
For i = 1 To MaxGeneration
For j = 1 To Popsize
ReDim InBird(i, j).v(VaryNum) As Double
ReDim InBird(i, j).x(VaryNum) As Double
Next j
Next i
For i = 1 To Popsize
ReDim IndividualBest(i).v(VaryNum) As Double
ReDim IndividualBest(i).x(VaryNum) As Double
Next i
'##########重定义数组变量##############
Randomize
'##########初始化微粒群##############
For i = 1 To Popsize
For j = 1 To VaryNum
InBird(1, i).x(j) = amin(j) + Rnd * (amax(j) - amin(j)) '初始化变量值
Next j
Next i
For i = 1 To Popsize
InBird(1, i).Value = FunctionObject(InBird(1, i).x(), 1) '通过调用函数计算目标函数值
IndividualBest(i) = InBird(1, i) '初始化时将初始值作为最优值
Next i
GlobalBest = IndividualBest(1)
For i = 1 To Popsize
If ObjectValue = 1 Then '目标函数取最大值为目的
If GlobalBest.Value < IndividualBest(i).Value Then
GlobalBest = IndividualBest(i) '寻求初始群体中的最优位置
End If
ElseIf ObjectValue = 0 Then '目标函数取最小值为目的
If GlobalBest.Value > IndividualBest(i).Value Then
GlobalBest = IndividualBest(i) '寻求初始群体中的最优位置
End If
End If
Next i
'##########初始化微粒群##############
For Generation = 2 To MaxGeneration
w = 1.2 * (1 - 0.9 * Generation / MaxGeneration) '0.4 - Generation * (0.8 - 0.6) / MaxGeneration '线性调整权重因子w的值
c1 = 0.5: c2 = 0.5 '比例参数
For i = 1 To Popsize
For j = 1 To VaryNum
Recalculation: '重新计算
Rnd1 = Rnd: Rnd2 = Rnd '产生随机数
'产生下一代微粒移动的距离
InBird(Generation, i).v(j) = w * InBird(Generation - 1, i).v(j) + c1 * Rnd1 * (IndividualBest(i).x(j) - InBird(Generation - 1, i).x(j)) + c2 * Rnd2 * (GlobalBest.x(j) - InBird(Generation - 1, i).x(j))
InBird(Generation, i).x(j) = InBird(Generation - 1, i).x(j) + InBird(Generation, i).v(j)
If InBird(Generation, i).x(j) > amax(j) Or InBird(Generation, i).x(j) < amin(j) Then
PNum = PNum + 1
If PNum <= 10 Then
GoTo Recalculation
Else
InBird(Generation, i).v(j) = Rnd
InBird(Generation, i).x(j) = amin(j) + (amax(j) - amin(j)) * InBird(Generation, i).v(j)
End If
End If
PNum = 0
Next j
Next i
For i = 1 To Popsize
InBird(Generation, i).Value = FunctionObject(InBird(Generation, i).x(), Generation)
Next i
For i = 1 To Popsize '寻求每一个微粒群个体运行至某一代时发现的最优位置
If ObjectValue = 1 Then '目标函数取最大值为目的
If IndividualBest(i).Value < InBird(Generation, i).Value Then
IndividualBest(i) = InBird(Generation, i)
End If
ElseIf ObjectValue = 0 Then '目标函数取最小值为目的
If IndividualBest(i).Value > InBird(Generation, i).Value Then
IndividualBest(i) = InBird(Generation, i)
End If
End If
Next i
For i = 1 To Popsize
If ObjectValue = 1 Then '目标函数取最大值为目的
If GlobalBest.Value < IndividualBest(i).Value Then
GlobalBest = IndividualBest(i) '寻求整个微粒群体运行至某一代时发现的最优位置
End If
ElseIf ObjectValue = 0 Then '目标函数取最小值为目的
If GlobalBest.Value > IndividualBest(i).Value Then
GlobalBest = IndividualBest(i) '寻求整个微粒群体运行至某一代时发现的最优位置
End If
End If
Next i
Debug.Print "The Current Generation is at"; Generation - 1; " ";
Debug.Print "GlobalBest.Value="; GlobalBest.Value
Next Generation
For i = 1 To VaryNum
Debug.Print "GlobalBest.x("; i; ")="; GlobalBest.x(i)
Next i
'Call SubValue(GlobalBest.x(), RealValue())
End Function
Public Sub SubValue(Vary() As Double, fValue() As Double)
Dim Dimension As Integer, x() As Double, i%
Dimension = UBound(Vary)
ReDim x(Dimension) As Double
For i = 1 To Dimension
x(i) = Vary(i)
Next i
'以下目标函数变时需要修改
Dim j%
For j = 1 To DataNum
fValue(j) = x(1) + x(2) * X0(j, 1) + x(3) * X0(j, 2) + x(4) * X0(j, 3) + x(5) * X0(j, 1) * X0(j, 2) + x(6) * X0(j, 1) * X0(j, 3) + x(7) * X0(j, 2) * X0(j, 3) + x(8) * X0(j, 1) ^ 2 + x(9) * X0(j, 2) ^ 2 + x(10) * X0(j, 3) ^ 2
Next j
'以上目标函数变时需要修改
For i = 1 To DataNum
Debug.Print "fValue("; i; ")="; fValue(i)
Next i
End Sub
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -