文化算法

基本概念：优化算法 | 详解文化算法（附MATLAB代码） - 知乎

不同于遗传算法只有种群进化空间，文化算法包含信念空间、种群空间两个进化空间，因此，文化算法是一种双层进化系统。

种群空间和信念空间是各自保存自己群体的两个相对独立的进化过程，并各自独立演化。但是，种群空间定期贡献精英个体给信念空间，信念空间不断进化自己的精英群体来影响种群空间群体。

信仰空间中的知识通过影响函数指导种群空间个体朝有利的方向进化。任何满足文化算法结构的算法都能融入种群空间中，通过相互融合，可以弥补传统算法的不足。

种群空间结构-适应度值

信念空间由两种知识组成：

1）每次迭代中发现最佳点的个体的确切位置。

2）保存和修正目标函数主变量的变化区间在可能发现最好解的范围内。

接收函数accept()的作用是将个体经验传递到信念空间，主要是在当前种群空间中选取优势个体，为信念空间的进化提供基础。接收函数一般为选取当前种群中最好的前nAccept个个体来更新信念空间。

更新函数：在接收到当前种群中最好的前nAccept个个体后，使用更新函数update()来更新信念空间。

影响函数：信念空间在形成群体经验后，通过influence()函数对群体空间中个体的行为规则进行修改，以使个体空间得到更高的进化效率。

主要代码

round 四舍五入为最近的小数或整数

Y = round(X) 将 X 的每个元素四舍五入为最近的整数。在对等情况下，即有元素的小数部分恰为 0.5 时，round 函数会偏离零四舍五入到具有更大幅值的整数。

Y = round(X,N) 四舍五入到 N 位数：

N > 0：舍入到小数点右侧的第 N 位数。

N = 0：四舍五入到最接近的整数。

N < 0：舍入到小数点左侧的第 N 位数。

Y = round(X,N,type) 指定四舍五入的类型。指定 'significant' 以四舍五入为 N 位有效数（从最左位数开始计数）。在此情况下，N 必须为正整数。

Y = round(t) 将 duration 数组 t 的每个元素四舍五入到最接近的秒数。

Y = round(t,unit) 将 t 的每个元素四舍五入到指定单位时间的最接近的数。

repmat 重复数组副本

B = repmat(A,n) 返回一个数组，该数组在其行维度和列维度包含 A 的 n 个副本。A 为矩阵时，B 大小为 size(A)*n。

B = repmat(A,r1,...,rN) 指定一个标量列表 r1,..,rN，这些标量用于描述 A 的副本在每个维度中如何排列。当 A 具有 N 维时，B 的大小为 size(A).*[r1...rN]。例如：repmat([1 2; 3 4],2,3) 返回一个 4×6 的矩阵。

B = repmat(A,r) 使用行向量 r 指定重复方案。例如，repmat(A,[2 3]) 与 repmat(A,2,3) 返回相同的结果。

inf 创建所有值均为inf的数组

X = Inf 返回正无穷大的标量表示。当运算结果太大以至于无法表示为浮点数时，如 1/0 或 log(0)，运算会返回 Inf。

对于双精度，Inf 表示大于 realmax 的数字。对于单精度，Inf 表示大于 realmax('single') 的数字。

X = Inf(n) 返回 Inf 值的 n×n 矩阵。

X = Inf(sz1,...,szN) 返回由 Inf 值组成的 sz1×...×szN 数组，其中 sz1,...,szN 指示每个维度的大小。例如：Inf(3,4) 返回一个 3×4 的矩阵。

X = Inf(sz) 返回 Inf 值的数组，其中大小向量 sz 定义 size(X)。例如：Inf([3 4]) 返回一个 3×4 的矩阵。

X = Inf(___,typename) 返回由数据类型为 typename（可以是 'single' 或 'double'）的 Inf 值组成的数组。

X = Inf(___,'like',p) 将返回一个由 Inf 值组成的数组，它具有与 p 相同的数据类型、稀疏度和复/实性。您可以指定 typename 或 'like'，但不能同时指定二者。

unifrnd 连续统一的随机数

r = unifrnd(a,b) 从连续均匀分布中生成一个随机数，其下端点为a，上端点为b。

r = unifrnd(a,b,sz1,...,szN) 生成一个均匀随机数的数组，其中sz1,...,szN表示每个维度的大小。

r = unifrnd(a,b,sz) 生成一个均匀随机数的数组，其中大小向量sz指定了size(r)。

semilogy 半对数图（y 轴有对数刻度）

semilogy(X,Y) 在 x 轴上使用线性刻度、在 y 轴上使用以 10 为底的对数刻度来绘制 x 和 y 坐标。

要绘制由线段连接的一组坐标，请将 X 和 Y 指定为相同长度的向量。

要在同一组坐标轴上绘制多组坐标，请将 X 或 Y 中的至少一个指定为矩阵。

semilogy(X,Y,LineSpec) 使用指定的线型、标记和颜色创建绘图。

semilogy(X1,Y1,...,Xn,Yn) 在同一组坐标轴上绘制多对 x 和 y 坐标。此语法可替代将坐标指定为矩阵的形式。

semilogy(X1,Y1,LineSpec1,...,Xn,Yn,LineSpecn) 可为每个 x-y 对组指定特定的线型、标记和颜色。您可以对某些 x-y 对组指定 LineSpec，而对其他对组省略它。例如，semilogy(X1,Y1,'o',X2,Y2) 对第一个 x-y 对组指定标记，但没有对第二个对组指定标记。

clc;
clear;
close all;

%% Problem Definition

CostFunction = @(x) Sphere(x);        % Cost Function

nVar = 10;          % Number of Decision Variables

VarSize = [1 nVar];   % Decision Variables Matrix Size

VarMin = -10;         % Decision Variables Lower Bound
VarMax = 10;         % Decision Variables Upper Bound

%% Cultural Algorithm Settings

MaxIt = 1000;         % Maximum Number of Iterations

nPop = 50;            % Population Size

pAccept = 0.35;                   % Acceptance Ratio接受的更新率
nAccept = round(pAccept*nPop);    % Number of Accepted Individuals

alpha = 0.3;

beta = 0.5;

%% Initialization

% Initialize Culture
Culture.Situational.Cost = inf;
Culture.Normative.Min = inf(VarSize);
Culture.Normative.Max = -inf(VarSize);
Culture.Normative.L = inf(VarSize);
Culture.Normative.U = inf(VarSize);

% Empty Individual Structure
empty_individual.Position = [];
empty_individual.Cost = [];

% Initialize Population Array
pop = repmat(empty_individual, nPop, 1);

% Generate Initial Solutions
for i = 1:nPop
    pop(i).Position = unifrnd(VarMin, VarMax, VarSize);
    pop(i).Cost = CostFunction(pop(i).Position);
end

% Sort Population
[~, SortOrder] = sort([pop.Cost]);
pop = pop(SortOrder);

% Adjust Culture using Selected Population
spop = pop(1:nAccept);
Culture = AdjustCulture(Culture, spop);

% Update Best Solution Ever Found
BestSol = Culture.Situational;

% Array to Hold Best Costs
BestCost = zeros(MaxIt, 1);

%% Cultural Algorithm Main Loop

for it = 1:MaxIt
    
    % Influnce of Culture
    for i = 1:nPop
        
        % % 1st Method (using only Normative component)
%         sigma = alpha*Culture.Normative.Size;
%         pop(i).Position = pop(i).Position+sigma.*randn(VarSize);
        
        % % 2nd Method (using only Situational component)
%         for j = 1:nVar
%            sigma = 0.1*(VarMax-VarMin);
%            dx = sigma*randn;
%            if pop(i).Position(j)<Culture.Situational.Position(j)
%                dx = abs(dx);
%            elseif pop(i).Position(j)>Culture.Situational.Position(j)
%                dx = -abs(dx);
%            end
%            pop(i).Position(j) = pop(i).Position(j)+dx;
%         end
        
        % % 3rd Method (using Normative and Situational components)
        for j = 1:nVar
          sigma = alpha*Culture.Normative.Size(j);
          dx = sigma*randn;
          if pop(i).Position(j)<Culture.Situational.Position(j)
              dx = abs(dx);
          elseif pop(i).Position(j)>Culture.Situational.Position(j)
              dx = -abs(dx);
          end
          pop(i).Position(j) = pop(i).Position(j)+dx;
        end        
        
        % % 4th Method (using Size and Range of Normative component)
%         for j = 1:nVar
%           sigma = alpha*Culture.Normative.Size(j);
%           dx = sigma*randn;
%           if pop(i).Position(j)<Culture.Normative.Min(j)
%               dx = abs(dx);
%           elseif pop(i).Position(j)>Culture.Normative.Max(j)
%               dx = -abs(dx);
%           else
%               dx = beta*dx;
%           end
%           pop(i).Position(j) = pop(i).Position(j)+dx;
%         end        
        
        pop(i).Cost = CostFunction(pop(i).Position);
        
    end
    
    % Sort Population
    [~, SortOrder] = sort([pop.Cost]);
    pop = pop(SortOrder);

    % Adjust Culture using Selected Population
    spop = pop(1:nAccept);
    Culture = AdjustCulture(Culture, spop);

    % Update Best Solution Ever Found
    BestSol = Culture.Situational;
    
    % Store Best Cost Ever Found
    BestCost(it) = BestSol.Cost;
    
    % Show Iteration Information
    disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);
    
end

%% Results

figure;
%plot(BestCost, 'LineWidth', 2);
semilogy(BestCost, 'LineWidth', 2);
xlabel('Iteration');
ylabel('Best Cost');
grid on;

Sphere

function z = Sphere(x)

    z = sum(x.^2);

end

AdjustCulture

numel - 数组元素的数目
此 MATLAB 函数返回数组 A 中的元素数目 n 等同于 prod(size(A))。

prod - 数组元素的乘积
此 MATLAB 函数返回 A 的数组元素的乘积。如果 A 是向量，则 prod(A) 返回元素的乘积。如果 A 为非空矩阵，则 prod(A) 将 A 的各列视为向量，并返回一个包含每列乘积的行向量。如果 A 为 0×0 空矩阵，prod(A) 返回 1。如果 A 为多维数组，则 prod(A) 沿第一个非单
一维度运算并返回乘积数组。此维度的大小将减少至 1，而所有其他维度的大小保持不变。

function Culture = AdjustCulture(Culture, spop)

    n = numel(spop);
    nVar = numel(spop(1).Position);
    
    for i = 1:n
        if spop(i).Cost<Culture.Situational.Cost
            Culture.Situational = spop(i);
        end
        
        for j = 1:nVar
            if spop(i).Position(j)<Culture.Normative.Min(j) ...
                    || spop(i).Cost<Culture.Normative.L(j)
                Culture.Normative.Min(j) = spop(i).Position(j);
                Culture.Normative.L(j) = spop(i).Cost;
            end
            if spop(i).Position(j)>Culture.Normative.Max(j) ...
                    || spop(i).Cost<Culture.Normative.U(j)
                Culture.Normative.Max(j) = spop(i).Position(j);
                Culture.Normative.U(j) = spop(i).Cost;
            end
        end
    end

    Culture.Normative.Size = Culture.Normative.Max-Culture.Normative.Min;
    
end