form1.frm

复合形法的不等式约束优化。其中的主程序是通用的

VERSION 5.00
Begin VB.Form Form1 
   AutoRedraw      =   -1  'True
   Caption         =   "Form1"
   ClientHeight    =   6840
   ClientLeft      =   60
   ClientTop       =   450
   ClientWidth     =   9000
   LinkTopic       =   "Form1"
   ScaleHeight     =   6840
   ScaleWidth      =   9000
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   615
      Left            =   3360
      TabIndex        =   0
      Top             =   6120
      Width           =   2295
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Option Explicit
    Dim Bmax As Double, Bmin As Double, G_ezai As Double
    Dim B_lg As Double, B_qb As Double, D1_qb As Double, D2_qb As Double
    Dim Beita As Double, Phi As Double
    Dim Ns As Double, Miu As Double, Sigma As Double
Private Sub Command1_Click()
    Dim i As Integer, N As Integer, M As Integer
    Dim alpha As Double, eps As Double, a() As Double, b() As Double, x() As Double, xx() As Double
    N = 3: M = 3    'N为变量个数；m为约束个数
    ReDim a(N - 1)   'a为常量约束条件中变量的下界，长度为N
    ReDim b(N - 1)   'b为a为常量约束条件中变量的上界，长度为N
    ReDim x(N)       '长度为n+1，前n个分量存放初始坐标，最后一个分量返回极小值
    ReDim xx((N + 1) * 2 * N - 1) '长度（n+1)*2n，最后一列的前n行值为最优点的坐标，最后一个值为最优点的函数值
   
    Bmax = 1.25            '最大工件宽度
    Bmin = 0.9             '最小工件宽度
    G_ezai = 8000            '夹具载荷
    Beita = 20           '夹取最大工件时钳臂与水平面的夹角
    Phi = 20             '夹取最大工件最小工件时钳臂的旋转角度
    Ns = 1.1              '预紧力的安全系数
    Miu = 0.35            '摩擦系数
    B_lg = 0.02             '连杆厚度
    B_qb = 0.03            '钳臂厚度
    Sigma = 100 * 10 ^ 6  '许用应力
    D1_qb = 0.08          '钳臂o点圆孔内径
    D2_qb = 0.06          '钳臂与连杆相连圆孔内径
    
    x(0) = 50: x(1) = 1.4: x(2) = 0.15 '初始值
    
    a(0) = 40: a(1) = (Bmax + Bmin) / 2: a(2) = D1_qb / 2
    b(0) = 90: b(1) = 2: b(2) = 0.5
    eps = 1# * 10 ^ -40: alpha = 1.3
    
    Randomize '初始化随机数

    i = Hecmp(N, M, a, b, alpha, eps, x, xx(), 500)   'i返回迭代次数
    Me.Print x(0)
    Me.Print x(1)
    Me.Print x(2)
    Me.Print x(3)
    Me.Print (i)

End Sub
'优化主函数
Public Function Hecmp(ByVal N As Integer, ByVal M As Integer, _
                 a() As Double, b() As Double, ByVal alpha As Double, _
                ByVal eps As Double, x() As Double, xx() As Double, ByVal K As Integer) As Integer
    
    Dim rn() As Double
    Dim r As Integer, G As Integer, i As Integer, j As Integer, it As Integer, kt As Integer, jt As Integer, kk As Integer
    Dim fj As Double, fr As Double, fg As Double, z As Double, rr As Double
    Dim C() As Double, D() As Double, W() As Double, xt() As Double, xf() As Double
    
    Dim mm As Double
    
    ReDim C(M - 1)
    ReDim D(M - 1)
    ReDim W(M - 1)
    ReDim xt(N - 1)
    ReDim xf(N - 1)
    
    '利用随机数产生可用的初始点
    Do
        Hecmpfs N, M, x(), C(), D(), W()
        r = 0: G = 0
        While ((r < N) And G = 0) '判断初始点x是否满足约束条件
            If ((C(r) <= W(r)) And (D(r) >= W(r))) Then
                r = r + 1
            Else
                G = 1
                For i = 0 To N - 1 '利用随机数产生初始点
                    x(i) = a(i) + (b(i) - a(i)) * Rnd
                Next i
            End If
        Wend
    Loop While G = 1
        
    
    
    rr = 0#
    'Randomize '初始化随机数
    For i = 0 To N - 1
        xx(i * N * 2) = x(i)
    Next i
    
    xx(N * N * 2) = Hecmpf(x(), N)
    
    For j = 1 To 2 * N - 1
        For i = 0 To N - 1
'            rr = 2053# * rr + 13849#
'            mm = rr / 65536#
'            rr = rr - mm * 65536#
'            xx(i * N * 2 + j) = a(i) + (b(i) - a(i)) * (rr / 65536#)
            'Randomize '初始化随机数
            xx(i * N * 2 + j) = a(i) + (b(i) - a(i)) * Rnd
            x(i) = xx(i * N * 2 + j)
        Next i
        
        it = 1
        
        While it = 1
            it = 0: r = 0: G = 0
            While ((r < N) And G = 0)
                If ((a(r) <= x(r)) And (b(r) >= x(r))) Then
                    r = r + 1
                Else
                    G = 1
                End If
            Wend
            
            If (G = 0) Then
                Hecmpfs N, M, x(), C(), D(), W()
                r = 0
                While ((r < M) And (G = 0))
                    If ((C(r) <= W(r)) And (D(r) >= W(r))) Then
                        r = r + 1
                    Else
                        G = 1
                    End If
                Wend
            End If
            
            If (G <> 0) Then
                For r = 0 To N - 1
                    z = 0#
                    For G = 0 To j - 1
                        z = z + xx(r * N * 2 + G) / (1# * j)
                    Next G
                    xx(r * N * 2 + j) = (xx(r * N * 2 + j) + z) / 2#
                    x(r) = xx(r * N * 2 + j)
                Next r
                it = 1
            Else
                xx(N * N * 2 + j) = Hecmpf(x, N)
            
            End If
        Wend
    Next j
    
    kk = 1: it = 1
    
    While (it = 1)
        it = 0
        fr = xx(N * N * 2): r = 0
        For i = 1 To 2 * N - 1
            If (xx(N * N * 2 + i)) > fr Then r = i: fr = xx(N * N * 2 + i)
        Next i
        G = 0: j = 0: fg = xx(N * N * 2)
        If r = 0 Then
            G = 1: j = 1: fg = xx(N * N * 2 + 1)
        End If
        For i = j + 1 To 2 * N - 1
            If (i <> r) Then
                If (xx(N * N * 2 + i) > fg) Then
                    G = i: fg = xx(N * N * 2 + i)
                End If
            End If
        Next i
        For i = 0 To N - 1
            xf(i) = 0#
            For j = 0 To 2 * N - 1
                If j <> r Then
                    xf(i) = xf(i) + xx(i * N * 2 + j) / (2# * N - 1#)
                End If
            Next j
            xt(i) = (1# + alpha) * xf(i) - alpha * xx(i * N * 2 + r)
        Next i
        jt = 1
        While jt = 1
            jt = 0
            z = Hecmpf(xt, N)
            While z > fg
                For i = 0 To N - 1
                    xt(i) = (xt(i) + xf(i)) / 2#
                Next i
                z = Hecmpf(xt, N)
            Wend
            j = 0
            For i = 0 To N - 1
                If (a(i) > xt(i)) Then
                    xt(i) = xt(i) + 0.0001: j = 1
                End If
                If (b(i) < xt(i)) Then
                    xt(i) = xt(i) - 0.0001: j = 1
                End If
            Next i
            If j <> 0 Then
                jt = 1
            Else
                Hecmpfs N, M, xt(), C(), D(), W()
                j = 0: kt = 1
                While ((kt = 1) And (j < M))
                    If ((C(j) <= W(j)) And (D(j) >= W(j))) Then
                        j = j + 1
                    Else
                        kt = 0
                    End If
                Wend
                If (j < M) Then
                    For i = 0 To N - 1
                        xt(i) = (xt(i) + xf(i)) / 2#
                    Next i
                    jt = 1
                End If
            End If
        Wend
        For i = 0 To N - 1
            xx(i * N * 2 + r) = xt(i)
        Next i
        xx(N * N * 2 + r) = z
        fr = 0#: fg = 0#
        For j = 0 To 2 * N - 1
            fj = xx(N * N * 2 + j)
            fr = fr + fj / (2# * N)
            fg = fg + fj * fj
        Next j
        fr = (fg - 2# * N * fr * fr) / (2# * N - 1#)
        If (fr >= eps) Then
            kk = kk + 1
            If (kk < K) Then it = 1
        End If
    Wend
    
    For i = 0 To N - 1
        x(i) = 0#
        For j = 0 To 2 * N - 1
            x(i) = x(i) + xx(i * N * 2 + j) / (2# * N)
        Next j
    Next i
    
    z = Hecmpf(x, N): x(N) = z
    
    Hecmpfs N, M, x(), C(), D(), W()
    
    Hecmp = kk
    
End Function
'目标函数
Private Function Hecmpf(x() As Double, N As Integer) As Double
    Dim y As Double, p As Double, q As Double, R1 As Double
    Dim Lab As Double, R2 As Double, H_lg As Double, Aa As Double, Bb As Double, Cc As Double
    H_lg = G_ezai * 9.8 / (8 * Cosdu(x(0)) * Sigma * B_lg)       '根据连杆受拉满应力计算连杆高度
    R2 = G_ezai * 9.8 / 16 / B_qb / Sigma / Cosdu(x(0)) + D2_qb / 2   '根据受拉计算与连杆相连处钳臂孔的外径
    'Aa = 8 * Sigma * B_qb * Cosdu(x(0))                                      '关于钳臂B点的外径的三次方程的系数
    'Bb = -1.5 * G_ezai * 9.8 * x(1) * Sindu(90 - x(0) + Beita)              '关于钳臂B点的外径的三次方程的系数
    'Cc = -Cosdu(x(0)) * Sigma * B_qb * D1_qb ^ 3                                     '关于钳臂B点的外径的三次方程的系数
    'p = Bb / Aa
    'q = Cc / Aa
    'R1 = Cuberoot(-q / 2 + Sqr(Abs((q / 2) ^ 2 + (p / 3) ^ 3))) + Cuberoot(-q / 2 - Sqr(Abs((q / 2) ^ 2 + (p / 3) ^ 3)))
    Lab = (x(1) * Cosdu(Beita) - (Bmax + Bmin) / 4) / Sindu(x(0)) '计算连杆AB的长度
    y = (Lab * 2 * B_lg * H_lg + (x(2) + R2) * x(1) * 2 * B_qb) '利用平均面积计算连杆和钳臂的体积
    Hecmpf = y                              '返回体积作为优化目标
End Function
'不等式约束
Private Sub Hecmpfs(N As Integer, M As Integer, x() As Double, C() As Double, D() As Double, W() As Double)
    'C不等式约束的下界
    'D不等式约束的上界
    'W不等式约束表达式
    Dim L1 As Double, n1 As Double, Lab As Double, Alpha2 As Double, Sinalpha2 As Double
    Dim L2 As Double, n2 As Double, M_qb As Double, W_qb As Double
    
    '条件0：夹取最大工件时加紧力校核
    L1 = x(1) * Sindu(90 - x(0) + Beita)  '轴o1到连杆AB的距离
    n1 = Tan(0.5 * DuToHudu(Phi)) * (1 + 4 * L1 / (Bmax - Bmin) / Cosdu(x(0))) '夹最大工件时的预紧力系数
    W(0) = n1 * Miu                              '夹取最大工件的预紧力约束条件方程
    C(0) = Ns                                     '该条件方程下界
    D(0) = 1 * 10 ^ 10                           '条件方程的上界
    
    '条件1：夹取最小工件时加紧力校核
    Sinalpha2 = Sindu(x(0)) * (x(1) * Cosdu(Beita + Phi) - (Bmax + Bmin) / 4) / (x(1) * Cosdu(Beita) - (Bmax + Bmin) / 4)
    If Sinalpha2 >= 1 Then
        Sinalpha2 = 0.9999
    End If
    If Sinalpha2 <= -1 Then
        Sinalpha2 = -0.9999
    End If
    Alpha2 = Atn(Sinalpha2 / Sqr(1 - Sinalpha2 * Sinalpha2))
    L2 = x(1) * Sindu(90 - HuduTodu(Alpha2) + Beita + Phi)
    n2 = Tan(0.5 * DuToHudu(Phi)) * (-1 + 4 * L2 / (Bmax - Bmin) / Cos(Alpha2))
    W(1) = n2 * Miu                              '夹取最大工件的预紧力约束条件方程
    C(1) = Ns                                     '该条件方程下界
    D(1) = 1 * 10 ^ 10                           '条件方程的上界
    
    '条件2：连杆受拉强度条件
    'W(2) = G_ezai * 9.8 / (8 * Cosdu(x(0)) * x(4) * B_lg)
    'C(2) = -1 * 10 ^ 10
    'D(2) = Sigma
    
    '条件3：钳臂o点受弯强度条件
    M_qb = 0.25 * G_ezai * 9.8 * x(1) * Sindu(90 - x(0) + Beita) / Cosdu(x(0))
    W_qb = B_qb * (8 * (x(2)) ^ 3 - D1_qb ^ 3) / 6 / x(2)
    W(2) = M_qb / W_qb
    C(2) = -1 * 10 ^ 10
    D(2) = Sigma
    
    '条件4：钳臂与连杆相连处抗拉强度条件
    'W(4) = G_ezai * 9.8 / 8 / B_qb / (2 * x(3) - D2_qb) / Cosdu(x(0))
    'C(4) = -1 * 10 ^ 10
    'D(4) = Sigma
        
End Sub


Private Function DuToHudu(du As Double) As Double
    DuToHudu = du / 180 * 3.1415926
End Function

Private Function HuduTodu(hudu As Double) As Double
    HuduTodu = hudu / 3.1415926 * 180
End Function
Private Function Sindu(du As Double) As Double
    Sindu = Sin(DuToHudu(du))
End Function

Private Function Cosdu(du As Double) As Double
    Cosdu = Cos(DuToHudu(du))
End Function
Private Function Cuberoot(xx As Double) As Double
    Dim croot As Double, root As Double
    croot = xx
    root = croot
    croot = (2 * root + xx / (root * root)) / 3
    While Abs(croot - root) > 0.000001
        root = croot
        croot = (2 * root + xx / (root * root)) / 3
    Wend
    Cuberoot = croot
End Function