将MTOA（多跟踪器优化算法）从MATLAB转成python的尝试

藜_

已于 2023-06-07 14:45:29 修改

阅读量192

点赞数

文章标签： matlab python 算法

于 2023-06-07 14:35:40 首次发布

本文链接：https://blog.csdn.net/qq_73347764/article/details/131086240

版权

总之就是写篇文章来记录下我被从没接触过的语言折磨的事情。如果有好心老哥愿意帮我改改的话……我会十分感激的。

学校开的人工智能入门的作业是找一个优化算法讲解，选的算法不能和别人重复，我选得有些晚了，原本是在mealpy里找算法的，但是算法库太小，最后在网上找到了Multi-tracker Optimization Algorithm(MTOA)——多跟踪器优化算法，奈何网上似乎只有MATLAB代码。我又不会这个，在GPT的帮助下，再加上我手动修修补补成现在这个样子了。

以下是在Springer上找到的MTOA代码：

原MATLAB文件“Creat_LTs”👇

function Out = Create_LTs(No_LTs,Rs,Dim)

R_p=Rs*(rand(1,No_LTs));
D_p=randn(Dim,No_LTs);
D_p=D_p./repmat(((sum(D_p.^2)).^(1/2)),Dim,1);
D_p(find(isnan(D_p))==1)=1;
Out=D_p.*repmat(R_p,Dim,1);

end

原MATLAB文件“Ev_Fcn”👇

function Out = Ev_Fcn(Points,Fcn_Name)

[m,n]=size(Points);

for i=1:n
o(1,i)=feval(Fcn_Name,Points(:,i));
end

Out=o;

end

原MATLAB文件“Gs”👇

function out = Gs( w,num,den,k,png,pdg )

GS_FOTF=k*fotf(den,pdg,num,png);
out=freqresp(GS_FOTF,w);

end

原MATLAB文件“MTOA”👇

function [OP_Cost,GOP]=MTOA(Fcn_Name,Par_Interval,No_GTs,No_LTs,RM,Rm,Max_Itr,Beta,Lambda,Theta,Graphic_on)

%                           initialization
%===============================================================
[Dim,m]=size(Par_Interval);
OP_Cost=zeros(1,Max_Itr);
LP=zeros(Dim,No_GTs);
LP_Cost=ones(1,No_GTs)*inf;

GTs=rand(Dim,No_GTs).*repmat(Par_Interval(:,2)-Par_Interval(:,1),1,No_GTs)+repmat(Par_Interval(:,1),1,No_GTs);
GTs_Cost=Ev_Fcn(GTs,Fcn_Name);
[Gts_Sorted,RKs]=sort(GTs_Cost);
GOP=GTs(:,RKs(1));
GOP_Cost=Gts_Sorted(1);
nop=No_GTs;
OP_Cost(1)=GOP_Cost;GOP
%===============================================================

for Itr=1:Max_Itr
    for i=1:nop
       %                Determin RSs and Search by LTs                  
       %--------------------------------------------------------
       Rf=((i-1)/(nop-1))*(RM-Rm)+Rm;
       Rd=norm(GOP-GTs(:,RKs(i)));
       Rs=Rf*(Rf>=Rd)+Rd*(Rd>Rf);
       
       LTs_C=Create_LTs(No_LTs,Rs,Dim);
       
       LTs=repmat(GTs(:,RKs(i)),1,No_LTs)+LTs_C;
       LTs=SS(LTs,Par_Interval);
       
       %----------------
       if Graphic_on==1
       subplot(2,2,1)
       hold off
       pause(0.000001);
       plot(LTs(1,:),LTs(2,:),'x');
       hold on
       ezplot(['(x-' num2str(GTs(1,RKs(i))) ')^2 + (y-' num2str(GTs(2,RKs(i))) ')^2 -' num2str(Rs^2)],[0 10],[0 10]);
       hold off
       xlim([Par_Interval(1,1) Par_Interval(1,2)]);
       ylim([Par_Interval(2,1) Par_Interval(2,2)]);
       pbaspect([1 1 1])
       title('Local Search')
       xlabel('x_1')
       ylabel('x_2')
       end
       %----------------
       
       LTs_Cost=Ev_Fcn(LTs,Fcn_Name);
       [L_min,L_inx]= min(LTs_Cost);
       if L_min<=LP_Cost(RKs(i))
        LP(:,RKs(i))=LTs(:,L_inx);
        LP_Cost(RKs(i))=L_min;
       end
       
       if L_min<=GOP_Cost
        GOP_Cost=L_min;
        GOP=LTs(:,L_inx);
       end
    
    end
    
    %                Search by GTs                  
    %--------------------------------------------------------
    for i=1:nop
       GTs(:,i)=New_GT(GTs(:,i),LP(:,i),GOP,Lambda,Theta,Beta);
       GTs(:,i)=SS(GTs(:,i),Par_Interval);
       GTs_Cost(i)=Ev_Fcn(GTs(:,i),Fcn_Name);
    end
    
    %                  Ranking
    %--------------------------------------------------------
    [Gts_Sorted,RKs]=sort(GTs_Cost);
    GOP_B=GTs(:,RKs(1));
    GOP_Cost_B=Gts_Sorted(1);
    if GOP_Cost_B<=GOP_Cost
        GOP_Cost=GOP_Cost_B;
        GOP=GOP_B;
    end
    OP_Cost(Itr+1)=GOP_Cost;
    
    %----------------
    if Graphic_on==1
    subplot(2,2,2)
    hold off
    pause(.000001)
    plot(GTs(1,:),GTs(2,:),'*')
    hold on
    plot(GOP(1,:),GOP(2,:),'X','color','red')
    xlim([Par_Interval(1,1) Par_Interval(1,2)]);
    ylim([Par_Interval(2,1) Par_Interval(2,2)]);
    hold off
    pbaspect([1 1 1]*3)
    title('Global Search')
    xlabel('x_1')
    ylabel('x_2')
    end
    %----------------
    %----------------
    if Graphic_on==1
    subplot(2,2,3)
    hold off
    pause(.000001)
    plot(OP_Cost(1:Itr+1))
    pbaspect([2 1 1])
    xlim([1 Max_Itr+1])
    title(['Cost=' num2str(GOP_Cost,'%4.10f')])
    xlabel('Iteration')
    ylabel('Cost')
    else
    hold off
    pause(.000001)
    plot(0:Itr,OP_Cost(1:Itr+1),'.','MarkerSize',15,'LineStyle','-','Color',[214 30 0]/255,'MarkerEdgeColor',[3 93 118]/255)
    pbaspect([2 1 1])
    title(['Itr=' num2str(Itr) ', Cost=' num2str(GOP_Cost,'%4.10f')]) 
    xlim([0 Max_Itr])
    xlabel('Iteration')
    ylabel('Cost')
    end
    %----------------
    
end


end

↑这个源文件不管怎么看都是主体部分↑

原MATLAB文件“New_GT”👇

function NGT = New_GT(GTs,LP,GOP,Lambda,Theta,Beta)


Dim=length(GTs);
G=Beta*(GOP-GTs)+(1-Beta)*(LP-GTs)+GTs;
d=(G-GTs);
D=norm(d);
IP=randn(Dim,1)*(D/10)+G;

%===============================
d2=(IP-GTs);

IP2=d*((d.'*d2)/(norm(d)^2))+GTs;
IP=IP-IP2;
Phi=Theta*rand;
Gamma=Lambda*rand;
IP=(IP/norm(IP))*tan(Phi)*D;
IP=IP+G;
V=(IP-GTs);
NGT=(V/norm(V))*D*Gamma+GTs;

if sum(isnan(NGT))
    NGT=GTs;
end

end

原MATLAB文件“SS”👇

function out= SS(Gs,Lim)

[m,n]=size(Gs);

for i=1:m
    for j=1:n
        Gs(i,j)=(Gs(i,j)<Lim(i,1))*Lim(i,1)+(Gs(i,j)>Lim(i,2))*Lim(i,2)+((Lim(i,1)<=Gs(i,j)) && (Gs(i,j)<=Lim(i,2)))*Gs(i,j);
    end
end

out=Gs;
end

原MATLAB文件“TestFcn”👇
function u = TestFcn(in)

x=in*2-1;
u=(x(1)^2+x(2)^2-cos(18*x(1))-cos(18*x(2)));

end
这是用作测试函数的公式。

原MATLAB文件“TestFcn”👇

clear
clc
close all
format long

%=======================Initializing MTOA Parameters=====================
Fcn_Name='TestFcn'; % Optimization Problrm Function
%=== Rastrigin Function (Function number 11 in the MTOA paper (Table. 2)) <----(Normalized Input)
%
%    x=xn*2-1;
%    u=(x(1)^2+x(2)^2-cos(18*x(1))-cos(18*x(2)));
%  
%    Min Cost = -2 @ xn=(0.5 , 0.5)

Min=zeros(1,2).'; 
Max=ones(1,2).'; 
Par_Interval=[Min Max]; %  Search Space Limitation 
No_GTs=20; %  Number of Global Tracker 
No_LTs=4; %  Number of Local Tracker
%---------- Equivalent Population = 100

RM=sqrt(2); % Maximum Search Radius
Rm=1e-4; % Minimum Search Radius
Max_Itr=100; % Maximum Iteration
Beta=.95; % Beta
Lambda=2; % Lambda
Theta=pi/8; % Theta

%==============================Graphic Option============================
TwoDGraphic_on=0; % (0  or  1) % If this option is enabled (TwoDGraphic_on=1) the Local and Global search process will be displayed.

%============================MTOA Function Run===========================
[GOP_Cost,GOP]=MTOA(Fcn_Name,Par_Interval,No_GTs,No_LTs,RM,Rm,Max_Itr,Beta,Lambda,Theta,TwoDGraphic_on);

%=========   ==================Solution==================================
disp('Solution=')
disp(GOP);

👇原文附带的测试函数生成图片

下面是目前的python代码，主要是GPT的功劳，我小改了一下，但还是有bug：


#调用库
import numpy as np
import control
import matplotlib.pyplot as plt


#Create_LTs（返回Out）
def create_LTs(No_LTs, Rs, Dim):    
    R_p = Rs * np.random.rand(1, No_LTs)    
    D_p = np.random.randn(Dim, No_LTs)    
    D_p = D_p / np.tile(np.sqrt(np.sum(D_p**2, axis=0)), (Dim, 1))    
    D_p[np.isnan(D_p)] = 1    
    Out = D_p * np.tile(R_p, (Dim, 1))    
    return Out

#Ev_Fcn（）
def Ev_Fcn(points, fcn_name):
    m, n = points.shape
    o = np.zeros((1, n))
    for i in range(n):
        o[0, i] = eval(fcn_name)(points[:, i])
    out = o
    return out

#Gs
def Gs(w, num, den, k, png, pdg):
    GS_FOTF = k*control.tf(num, den)*control.delay(pdg)/control.delay(png)
    out = control.freqresp(GS_FOTF, w)
    return out

#MTOA
def MTOA(Fcn_Name, Par_Interval, No_GTs, No_LTs, RM, Rm, Max_Itr, Beta, Lambda, Theta, Graphic_on):
    OP_Cost = None
    GOP = None
    # function body to be implemented

    #initialization
    Dim, m = Par_Interval.shape
    OP_Cost = np.zeros((1, Max_Itr))
    LP = np.zeros((Dim, No_GTs))
    LP_Cost = np.ones((1, No_GTs)) * np.inf

    GTs = np.random.rand(Dim, No_GTs) * np.tile(Par_Interval[:, 1] - Par_Interval[:, 0], (1, No_GTs)) + np.tile(Par_Interval[:, 0], (1,No_GTs))
    
    #bug
    
    
    GTs_Cost = Ev_Fcn(GTs, Fcn_Name)
    Gts_Sorted, RKs = np.sort(GTs_Cost), np.argsort(GTs_Cost)
    GOP = GTs[:, RKs[0]]
    GOP_Cost = Gts_Sorted[0]
    nop = No_GTs
    OP_Cost[0] = GOP_Cost
    '''
    Please note that the Ev_Fcn() function is not defined in the code snippet provided, so it needs to be implemented or imported from an external module.
    '''
    
    for Itr in range(1, Max_Itr+1):
        # loop over the range from 1 to nop
        for i in range(1, nop+1):
            # Determin RSs and Search by LTs
            Rf = ((i-1)/(nop-1))*(RM-Rm)+Rm
            Rd = np.linalg.norm(GOP-GTs[:,RKs[i]])
            Rs = Rf*(Rf>=Rd)+Rd*(Rd>Rf)
            
            LTs_C = create_LTs(No_LTs,Rs,Dim)  # create local search points
            
            LTs = np.tile(GTs[:,RKs[i]], (1, No_LTs)) + LTs_C
            LTs = np.SS(LTs, Par_Interval)  # do some operation
            
            # plot graph if Graphic_on is True
            if Graphic_on==1:
                plt.subplot(2,2,1)
                plt.hold(False)
                plt.pause(0.000001)
                plt.plot(LTs[0,:], LTs[1,:], 'x')
                plt.hold(True)
                plt.plot(['(x-',str(GTs(1,RKs(i))),')**2 + (y-',str(GTs(2,RKs(i))),')**2 -',str(Rs**2)],[0,10],[0,10])
                plt.hold(False)
                plt.xlim([Par_Interval[0][0], Par_Interval[0][1]])
                plt.ylim([Par_Interval[1][0], Par_Interval[1][1]])
                plt.pbaspect([1, 1, 1])
                plt.title('Local Search')
                plt.xlabel('x_1')
                plt.ylabel('x_2')
            
            LTs_Cost = Ev_Fcn(LTs, Fcn_Name)  # evaluate cost
            [L_min, L_inx] = np.min(LTs_Cost)
            
            # update cost and position of LP and GOP, if applicable
            if L_min <= LP_Cost[RKs[i]]:
                LP[:, RKs[i]] = LTs[:, L_inx]
                LP_Cost[RKs[i]] = L_min
            if L_min <= GOP_Cost:
                GOP_Cost = L_min
                GOP = LTs[:, L_inx]
        # Search by GTs
        for i in range(nop):
            GTs[:,i]=New_GT(GTs[:,i], LP[:,i], GOP, Lambda, Theta, Beta)
            GTs[:,i]=SS(GTs[:,i], Par_Interval)
            GTs_Cost[i]=Ev_Fcn(GTs[:,i], Fcn_Name)

        # Ranking
        Gts_Sorted, RKs = np.sort(GTs_Cost)
        GOP_B = GTs[:, RKs[0]]
        GOP_Cost_B = Gts_Sorted[0]
        if GOP_Cost_B <= GOP_Cost:
            GOP_Cost = GOP_Cost_B
            GOP = GOP_B
        OP_Cost[Itr+1] = GOP_Cost

        if Graphic_on == 1:
            plt.subplot(2, 2, 2)
            plt.hold(False)
            plt.pause(.000001)
            plt.plot(GTs[0,:], GTs[1,:], '*')
            plt.hold(True)
            plt.plot(GOP[0,:], GOP[1,:], 'X', color='red')
            plt.xlim([Par_Interval[0,0], Par_Interval[0,1]])
            plt.ylim([Par_Interval[1,0], Par_Interval[1,1]])
            plt.hold(False)
            plt.pbaspect([1, 1, 1] * 3)
            plt.title('Global Search')
            plt.xlabel('x_1')
            plt.ylabel('x_2')

        if Graphic_on == 1:
            plt.subplot(2,2,3)
            plt.hold(False)
            plt.pause(0.000001)
            plt.plot(OP_Cost[0:Itr+1])
            plt.gca().set_aspect([2,1,1])
            plt.xlim([1,Max_Itr+1])
            plt.title('Cost=' + '{:4.10f}'.format(GOP_Cost))
            plt.xlabel('Iteration')
            plt.ylabel('Cost')
        else:
            plt.hold(False)
            plt.pause(0.000001)
            plt.plot(range(Itr+1), OP_Cost[0:Itr+1], linestyle='-', marker='.', markersize=15, color=[214, 30, 0]/255, markeredgecolor=[3, 93, 118]/255)
            plt.gca().set_aspect([2,1,1])
            plt.title('Itr=' + str(Itr) + ', Cost=' + '{:4.10f}'.format(GOP_Cost))
            plt.xlim([0, Max_Itr])
            plt.xlabel('Iteration')
            plt.ylabel('Cost')
    return OP_Cost, GOP

#New_GT
def New_GT(GTs,LP,GOP,Lambda,Theta,Beta):
    Dim = len(GTs)
    G = Beta * (GOP - GTs) + (1 - Beta) * (LP - GTs) + GTs
    d = (G - GTs)
    D = np.linalg.norm(d)
    IP = np.random.randn(Dim, 1) * (D / 10) + G

    # ===============================
    d2 = (IP - GTs)
    IP2 = d * ((d.T @ d2) / (np.linalg.norm(d) ** 2)) + GTs
    IP = IP - IP2
    Phi = Theta * np.random.rand()
    Gamma = Lambda * np.random.rand()
    IP = (IP / np.linalg.norm(IP)) * np.tan(Phi) * D
    IP = IP + G
    V = (IP - GTs)
    NGT = (V / np.linalg.norm(V)) * D * Gamma + GTs

    if np.isnan(NGT).any():
        NGT = GTs

    return NGT


#SS
def SS(Gs, Lim):
    m, n = Gs.shape
    for i in range(m):
        for j in range(n):
            Gs[i, j] = (Gs[i, j] < Lim[i, 0]) * Lim[i, 0] \
                        + (Gs[i, j] > Lim[i, 1]) * Lim[i, 1] \
                        + ((Lim[i, 0] <= Gs[i, j]) and (Gs[i, j] <= Lim[i, 1])) * Gs[i, j]
    return Gs

#TestFcn
def TestFcn(in_):
    x = in_ * 2 - 1
    u = x[0]**2 + x[1]**2 - np.cos(18*x[0]) - np.cos(18*x[1])
    return u



#最后一段
# Define test function
def test_fcn(x):
    x = x*2 - 1
    u = x[0]**2 + x[1]**2 - np.cos(18*x[0]) - np.cos(18*x[1])
    return -u

# Set MTOA parameters
min_val = np.zeros((2, 1))
max_val = np.ones((2, 1))
par_interval = np.concatenate((min_val, max_val), axis=1)
no_gts = 20
no_lts = 4
rm = np.sqrt(2)
r_m = 1e-4
max_itr = 100
beta = 0.95
lmbda = 2
theta = np.pi / 8
two_d_graphic_on = 0

# Run MTOA function
gop_cost, gop = MTOA(test_fcn, par_interval, no_gts, no_lts, rm, r_m, max_itr, beta, lmbda, theta, two_d_graphic_on)

# Display solution
print("Solution=")
print(gop)

在spyder里运行了一下（当然没成功），报错是这样的：

runfile('C:/Users/GI/untitled0.py', wdir='C:/Users/GI')
Traceback (most recent call last):

File "C:\Users\GI\untitled0.py", line 214, in <module>
gop_cost, gop = MTOA(test_fcn, par_interval, no_gts, no_lts, rm, r_m, max_itr, beta, lmbda, theta, two_d_graphic_on)

File "C:\Users\GI\untitled0.py", line 44, in MTOA
GTs = np.random.rand(Dim, No_GTs) * np.tile(Par_Interval[:, 1] - Par_Interval[:, 0], (1, No_GTs)) + np.tile(Par_Interval[:, 0], (1,No_GTs))

ValueError: operands could not be broadcast together with shapes (2,20) (1,40)

总之应该是MTOA函数里面的“GTs = np.random.rand(Dim, No_GTs) * np.tile(Par_Interval[:, 1] - Par_Interval[:, 0], (1, No_GTs)) + np.tile(Par_Interval[:, 0], (1,No_GTs))”这一行出了问题。

以上。