混沌优化算法（COA）

内容

1. 引言
2. 算法的实现

引言

在现代计算问题中，尤其是金融市场与算法交易场景下，高效的优化方法起着至关重要的作用。在众多全局优化算法中，受自然现象启发而设计的算法占据着独特地位。混沌优化算法（COA）便是这类创新方案之一，它将混沌理论与全局优化方法融为一体。

本文提出一种改进型混沌优化算法，该算法充分利用了混沌系统的基本特性：遍历性、对初始条件的敏感性以及类随机行为。算法通过多种混沌映射高效遍历搜索空间，尽可能避免陷入局部最优 —— 这是优化算法中常见的难点问题。

该算法的显著特点，是将混沌搜索与加权梯度法相结合，并引入自适应机制，实现搜索参数的动态调整。通过采用逻辑映射、正弦映射、帐篷映射等多种混沌映射，算法在复杂多峰函数上表现出良好的抗停滞能力，且能够有效找到全局最优解。我们将照例先完整解析算法原理，再通过常用标准测试函数验证其性能。 

算法的实现

该算法主要包含三个阶段：

• 第一载波混沌搜索，初始化混沌变量，通过逻辑映射生成连续的混沌变量值，再将混沌变量映射到待优化变量的取值区间内，并保存搜索得到的最优解。
• 沿加权梯度方向搜索；
• 第二载波混沌搜索，在最优解附近开展局部搜索，并采用分形方式设置搜索步长。 

在下图中，我将直观展示混沌优化算法的核心思想：把混沌作为一种实用的优化工具，而非随机过程。全局最优解（中心位置的亮黄色光晕）是我们要寻找的目标。蓝色发光粒子代表搜索智能体，它们沿着混沌轨迹移动；这些轨迹以发光线条呈现，体现出运动的非线性特征。

混沌行为关键特性展示：确定性（轨迹平滑连贯，并非随机杂乱）；遍历性（粒子能够探索整个搜索空间）；对初始条件敏感（不同粒子沿不同轨迹移动）。同时，搜索动态过程以不同强度的光晕表现，反映出不同区域的搜索 “能量”。最优解周围的同心圆代表吸引域，模糊渐变与梯度效果则体现出搜索空间的连续性。算法的主要阶段：

• 远离中心的大范围搜索（外围粒子），
• 逐步向优质区域逼近（中等轨迹），
• 最优解附近的局部搜索（靠近中心的粒子）。

这张效果图可以看作是混沌优化的一幅 “画像”：混沌并非无序混乱，而是一种可控的解空间探索过程。

图例 1. 混沌优化算法可视化

下图所示的算法流程图反映了三个主要阶段：

1. 第一载波搜索（蓝色模块）：利用混沌映射进行全局搜索，并将混沌变量映射至搜索空间；
2. 加权梯度搜索（橙色模块）：采用加权梯度法，通过权重系数处理约束条件；
3. 第二载波搜索（紫色模块）：在最优解附近进行局部搜索，并对参数 α 进行自适应调整。           

 图例 2. 混沌优化算法运行原理图

该图展示了混沌优化算法（CHAOS/COA）的基础三阶段结构。在经过多次尝试后，我最终选定了当前这个更先进、更实用的版本：不仅扩展了智能体结构，还新增了移动计数器、停滞计数器，以及为每个维度单独设置的混沌变量。丰富混沌映射的种类也十分关键；原作者仅使用了逻辑映射，而我的版本额外加入了正弦映射和帐篷映射。本次实现还包含了惩罚参数的自适应调整机制，以及能够根据优化成效和惯性进行动态调节的搜索参数机制。此外，还加入了基于拉丁超立方抽样的更复杂初始化方式。

我们开始编写算法伪代码：

算法初始化

1. 设置算法参数：
   • 种群规模（popSize）
   • 第一阶段迭代次数（S1）
   • 第二阶段迭代次数（S2）
   • 惩罚参数（sigma）
   • 调整系数（t3）
   • 权重系数小常数（eps）
   • 惯性参数（inertia）
   • 社会影响系数（socialFactor）
   • 变异概率（mutationRate）
2. 智能体初始化：
   • 对种群中的每个智能体：
     • 使用不同初始值初始化混沌变量（gamma）
     • 重置速度向量（velocity）
     • 将停滞计数器置为 0
3. 搜索参数（alpha）初始化：
   • 根据搜索空间大小自适应调整参数
   • alpha[c] = 0.1 * (rangeMax[c] - rangeMin[c]) / sqrt(coords)

第一阶段：基于混沌分布的初始种群构建

1. 使用拉丁超立方抽样生成初始种群：
   • 为每个维度生成数值
   • 打乱数值以保证均匀分布
   • 将数值映射到变量取值区间
2. 采用多种智能体初始化策略：
   • 前四分之一智能体：均匀分布
   • 第二个四分之一：围绕多个点聚类分布
   • 第三个四分之一：带边界偏移的随机位置
   • 最后四分之一：使用不同映射的混沌位置

第二阶段：第一载波混沌搜索

1. 对每一次迭代直至 S1：
   • 对种群中的每个智能体：
     • 若触发变异概率，则执行变异操作
     • 否则，对每个坐标：
       • 通过选定映射（逻辑映射、正弦映射或帐篷映射）更新混沌变量
       • 根据优化阶段确定搜索策略（全局搜索或局部搜索）
       • 全局搜索：使用混沌值确定新位置
       • 局部搜索：结合最优解吸引与混沌扰动
       • 结合惯性更新速度
       • 施加速度与位置约束
       • 检查约束越界并执行修正
   • 评估与更新：
     • 计算每个智能体的惩罚函数值
     • 更新个体最优解与全局最优解
     • 动态更新惩罚参数
     • 根据优化效果自适应调整搜索参数（alpha）

第三阶段：第二载波混沌搜索

1. 对从 S1 到 S1+S2 的每一次迭代：
   • 检查算法收敛性
   • 对种群中的每个智能体：
     • 若检测到收敛，则对部分智能体执行随机变异
     • 否则，对每个坐标：
       • 更新混沌变量
       • 自适应缩小搜索半径
       • 优先以最优解为基点选择基准位置
       • 添加混沌扰动与莱维噪声以实现随机长跳跃
       • 结合惯性更新速度与位置
       • 施加位置约束
   • 评估与更新：
     • 计算每个智能体的惩罚函数值
     • 更新个体最优解与全局最优解
     • 更新全局最优解历史以判断收敛性
     • 必要时重置停滞智能体

辅助函数

1. 混沌映射：
   • 逻辑映射 LogisticMap(x)：x [n+1] = r・x [n]・(1−x [n])，附带合法性校验
   • 正弦映射 SineMap(x)：x [n+1] = sin (π・x [n])，附带归一化处理
   • 帐篷映射 TentMap(x)：x [n+1] = μ・min (x [n], 1−x [n])，附带合法性校验
   • 选择混沌映射 SelectChaosMap(x, type)：根据类型选择对应映射
2. 约束与惩罚处理：
   • CalculateConstraintValue(agent, coord)：计算约束违反值
   • CalculateConstraintGradient(agent, coord)：计算约束梯度
   • CalculateWeightedGradient(agent, coord)：计算加权梯度
   • CalculatePenaltyFunction(agent)：计算惩罚函数
3. 停滞与收敛处理：
   • IsConverged()：基于最优解历史判断收敛
   • ResetStagnatingAgents()：重置停滞智能体
   • ApplyMutation(agent)：应用多种变异方式
   • UpdateSigma()：动态更新惩罚参数
   • UpdateBestHistory(newBest)：更新最优值历史

接下来我们开始讲解 COA（CHAOS）算法的具体实现。本文将详解所有核心实现方法，下一篇文章将直接进入算法测试与性能结果分析。我们来实现 S_COA_Agent 结构体及其字段：

• gamma []：一组伪随机类型的混沌变量，用于为智能体行为引入随机性与多样性，
• velocity []：速度数组，使智能体在搜索空间内更动态地移动，并考虑惯性作用，
• stagnationCounter：停滞计数器，当智能体解未出现优化时递增，用于实现停滞时的策略重启机制。

Init () 方法负责创建并设置数组初始值。gamma [] 采用 0.1 到 0.9 之间的均匀分布，为混沌变量的初始条件引入多样性。velocity [] 速度初始化为 0，停滞计数器置为 0。

//——————————————————————————————————————————————————————————————————————————————
// Improved agent structure with additional fields
struct S_COA_Agent
{
    double gamma    [];       // chaotic variables
    double velocity [];       // movement velocity (to increase inertia)
    int    stagnationCounter; // stagnation counter

    void Init (int coords)
    {
      ArrayResize (gamma,    coords);
      ArrayResize (velocity, coords);

      // Uniform distribution of gamma values
      for (int i = 0; i < coords; i++)
      {
        // Use different initial values for better variety
        gamma [i] = 0.1 + 0.8 * (i % coords) / (double)MathMax (1, coords - 1);

        // Initialize velocity with zeros
        velocity [i] = 0.0;
      }

      stagnationCounter = 0;
    }
};
//——————————————————————————————————————————————————————————————————————————————

C_AO_COA_chaos 类继承自基类 C_AO，是混沌优化算法（COA/CHAOS）的具体实现类。该类包含算法运行所需的方法与参数，同时基于混沌搜索理念，提供了用于控制智能体行为的扩展功能。类组成部分：

• SetParams() —— 设置算法参数的方法，
• Init() —— 初始化方法，用于接收算法运行所需的区间范围与参数，
• Moving() —— 负责智能体在解空间中移动的方法，
• Revision() —— 修正智能体位置的方法。

算法参数：

• S1、S2 —— 算法两个阶段的迭代次数。
• sigma、t3、eps、inertia、socialFactor、mutationRate —— 影响智能体行为与整体算法运行的参数。
• alpha[] —— 搜索过程使用的参数数组。
• agent[] —— 构成算法种群的智能体数组。

受保护方法：

• 梯度计算、约束值计算、惩罚函数计算，以及解可行性校验（IsFeasible）方法，
• 各类混沌映射方法（LogisticMap、SineMap、TentMap、SelectChaosMap）。

方法私有字段：

• 存储当前迭代周期信息的变量（epochs、epochNow）与动态惩罚值（currentSigma），
• globalBestHistory [] —— 用于存储多轮迭代中全局最优值的数组，
• 用于追踪最优值数组位置的历史索引（historyIndex），
• 管理智能体种群的方法（InitialPopulation）、执行不同搜索阶段的方法（FirstCarrierWaveSearch、SecondCarrierWaveSearch）、智能体变异方法（ApplyMutation）、惩罚更新方法（UpdateSigma）、收敛性检测方法（IsConverged），以及停滞智能体重置方法（ResetStagnatingAgents）。

综上，C_AO_COA_chaos 类是优化系统中的复杂核心组件，通过智能体实现最优解搜索。该类集成了控制算法内智能体所需的全部参数、方法与逻辑，同时包含确定性策略与混沌策略。 

//——————————————————————————————————————————————————————————————————————————————
class C_AO_COA_chaos : public C_AO
{
  public: //--------------------------------------------------------------------
  ~C_AO_COA_chaos () { }
  C_AO_COA_chaos ()
  {
    ao_name = "COA(CHAOS)";
    ao_desc = "Chaos Optimization Algorithm";
    ao_link = "https://www.mql5.com/en/articles/16729";

    // Internal parameters (not externally configurable)
    inertia      = 0.7;
    socialFactor = 1.5;
    mutationRate = 0.05;

    // Default parameters
    popSize = 50;
    S1      = 30;
    S2      = 20;
    sigma   = 2.0;
    t3      = 1.2;
    eps     = 0.0001;

    // Initialize the parameter array for the C_AO interface
    ArrayResize (params, 6);

    params [0].name = "popSize"; params [0].val = popSize;
    params [1].name = "S1";      params [1].val = S1;
    params [2].name = "S2";      params [2].val = S2;
    params [3].name = "sigma";   params [3].val = sigma;
    params [4].name = "t3";      params [4].val = t3;
    params [5].name = "eps";     params [5].val = eps;
  }

  void SetParams ()
  {
    // Update internal parameters from the params array
    popSize = (int)params [0].val;
    S1      = (int)params [1].val;
    S2      = (int)params [2].val;
    sigma   = params      [3].val;
    t3      = params      [4].val;
    eps     = params      [5].val;
  }

  bool Init (const double &rangeMinP  [], // minimum search range
             const double &rangeMaxP  [], // maximum search range
             const double &rangeStepP [], // search step
             const int     epochsP = 0);  // number of epochs

  void Moving   ();
  void Revision ();

  //----------------------------------------------------------------------------
  // External algorithm parameters
  int    S1;             // first phase iterations
  int    S2;             // second phase iterations
  double sigma;          // penalty parameter
  double t3;             // alpha correction ratio
  double eps;            // small number for weighting ratios

  // Internal algorithm parameters
  double inertia;        // inertia parameter for movement (internal)
  double socialFactor;   // social influence parameter (internal) 
  double mutationRate;   // mutation probability (internal) 

  S_COA_Agent agent [];  // array of agents

  private: //-------------------------------------------------------------------
  int    epochNow;
  double currentSigma;           // Dynamic penalty parameter
  double alpha             [];   // search parameters
  double globalBestHistory [10]; // History of global best solution values
  int    historyIndex;

  // Auxiliary methods
  double CalculateWeightedGradient   (int agentIdx, int coordIdx);
  double CalculateConstraintValue    (int agentIdx, int coordIdx);
  double CalculateConstraintGradient (int agentIdx, int coordIdx);
  double CalculatePenaltyFunction    (int agentIdx);

  // Method for checking the solution feasibility
  bool IsFeasible              (int agentIdx);

  // Chaotic maps
  double LogisticMap           (double x);
  double SineMap               (double x);
  double TentMap               (double x);
  double SelectChaosMap        (double x, int type);

  void InitialPopulation       ();
  void FirstCarrierWaveSearch  ();
  void SecondCarrierWaveSearch ();
  void ApplyMutation           (int agentIdx);
  void UpdateSigma             ();
  void UpdateBestHistory       (double newBest);
  bool IsConverged             ();
  void ResetStagnatingAgents   ();
};
//——————————————————————————————————————————————————————————————————————————————

C_AO_COA 类的 Init 方法负责算法的初始配置与运行前准备工作。它会执行多项关键任务：首先通过 StandardInit() 方法完成基础初始化，该方法用于设置搜索区间、步长及其他参数。若基础初始化执行失败，该方法会直接终止并返回错误。

随后，方法会设置与迭代总次数、当前迭代轮次（epochNow）以及惩罚系数（currentSigma）相关的参数。同时初始化最优解历史记录，用于追踪算法优化进程。方法会校验搜索区间的最小值数组与最大值数组长度是否一致。若数组长度不匹配或未指定区间，初始化流程将立即中断。

接下来，方法会初始化用于存储智能体、alpha 系数以及已找到最优解的数组。每个智能体都会根据不同策略生成初始位置：

• 一部分智能体在整个搜索区间内均匀分布，
• 另一部分采用聚类分布—— 在区间内的若干中心点附近分布，
• 其他智能体结合边界约束随机生成位置，
• 剩余智能体通过混沌映射函数生成初始解。

C_AO_COA_chaos 类的 Init 方法用于设置搜索最优解所需的初始参数与数组。该流程包括：校验输入数据的合法性、配置搜索区间、采用多种初始位置策略初始化智能体数组，以及设置已找到最优解等全局变量。方法执行期间，会创建后续迭代优化所需的数据结构，并设置用于调控智能体行为与整体算法运行的参数。

//——————————————————————————————————————————————————————————————————————————————
bool C_AO_COA_chaos::Init (const double &rangeMinP  [],
                     const double &rangeMaxP  [],
                     const double &rangeStepP [],
                     const int     epochsP = 0)
{
  if (!StandardInit (rangeMinP, rangeMaxP, rangeStepP)) return false;

  //----------------------------------------------------------------------------
  epochNow     = 0;
  currentSigma = sigma;
  historyIndex = 0;

  // Initialize the history of best values
  for (int i = 0; i < 10; i++) globalBestHistory [i] = -DBL_MAX;

  // Check and initialize the main arrays
  int arraySize = ArraySize (rangeMinP);
  if (arraySize <= 0 || arraySize != ArraySize (rangeMaxP) || arraySize != ArraySize (rangeStepP))
  {
    return false;
  }

  ArrayResize (agent, popSize);
  ArrayResize (alpha, coords);

  // Adaptive alpha initialization depending on the search range
  for (int c = 0; c < coords; c++)
  {
    // alpha depends on the size of the search space
    double range = rangeMax [c] - rangeMin [c];
    alpha [c] = 0.1 * range / MathSqrt (MathMax (1.0, (double)coords));
  }

  // Initialize of agents with various strategies
  for (int i = 0; i < popSize; i++)
  {
    agent [i].Init (coords);

    for (int c = 0; c < coords; c++)
    {
      double position;

      // Different initialization strategies
      if (i < popSize / 4)
      {
        // Uniform distribution in space
        position = rangeMin [c] + (i * (rangeMax [c] - rangeMin [c])) / MathMax (1, popSize / 4);
      }
      else
        if (i < popSize / 2)
        {
          // Clustering around multiple points
          int cluster = (i - popSize / 4) % 3;
          double clusterCenter = rangeMin [c] + (cluster + 1) * (rangeMax [c] - rangeMin [c]) / 4.0;
          position = clusterCenter + u.RNDfromCI (-0.1, 0.1) * (rangeMax [c] - rangeMin [c]);
        }
        else
          if (i < 3 * popSize / 4)
          {
            // Random positions with an offset towards the boundaries
            double r = u.RNDprobab ();
            if (r < 0.5) position = rangeMin [c] + 0.2 * r * (rangeMax [c] - rangeMin [c]);
            else position = rangeMax [c] - 0.2 * (1.0 - r) * (rangeMax [c] - rangeMin [c]);
          }
          else
          {
            // Chaotic positions using different maps
            int mapType = i % 3;
            double chaosValue = SelectChaosMap (agent [i].gamma [c], mapType);
            position = rangeMin [c] + chaosValue * (rangeMax [c] - rangeMin [c]);
          }

      a [i].cB [c] = u.SeInDiSp (position, rangeMin [c], rangeMax [c], rangeStep [c]);
    }
  }

  return true;
}
//——————————————————————————————————————————————————————————————————————————————

LogisticMap 方法实现了用于生成混沌序列的逻辑映射。该函数在算法中用于为解搜索过程引入随机性与多样性。该方法的核心思路是：基于当前状态值，通过逻辑映射公式计算新的状态值，并对公式参数进行微小调整，以强化混沌特性。

计算前，方法会校验输入值的有效性与区间合规性；若不符合要求，则用指定区间内的随机数替换。计算出新值后，同样会校验其是否在可接受范围内，必要时用随机数替换，以保证函数运行稳定。最终，方法通过内部逻辑确保生成的下一个状态值始终保持在允许的范围内。

//——————————————————————————————————————————————————————————————————————————————
// Improved chaotic maps
double C_AO_COA_chaos::LogisticMap (double x)
{
  // Protection against incorrect inputs
  if (x < 0.0 || x > 1.0 || MathIsValidNumber (x) == false)
  {
    x = 0.2 + 0.6 * u.RNDprobab ();
  }

  // x(n+1) = r*x(n)*(1-x(n))
  double r = 3.9 + 0.1 * u.RNDprobab (); // Slightly randomized parameter to avoid loops
  double result = r * x * (1.0 - x);

  // Additional validation
  if (result < 0.0 || result > 1.0 || MathIsValidNumber (result) == false)
  {
    result = 0.2 + 0.6 * u.RNDprobab ();
  }

  return result;
}
//——————————————————————————————————————————————————————————————————————————————

SineMap 方法实现了基于正弦函数的混沌映射。该方法接收当前状态值，先校验其有效性；如果数值不合法或超出 [0,1] 范围，就用该范围内的随机值进行替换。随后通过正弦函数计算新数值，并将其归一化回到 [0,1] 区间内，再执行一次额外校验。

如果最终数值越界或无效，会再次替换为 [0.2, 0.8] 区间内的随机数。最终，该方法返回基于当前状态、通过混沌映射生成的新状态值。

//——————————————————————————————————————————————————————————————————————————————
double C_AO_COA_chaos::SineMap (double x)
{
  // Protection against incorrect inputs
  if (x < 0.0 || x > 1.0 || MathIsValidNumber (x) == false)
  {
    x = 0.2 + 0.6 * u.RNDprobab ();
  }

  // x(n+1) = sin(π*x(n))
  double result = MathSin (M_PI * x);

  // Normalize the result to the range [0, 1]
  result = (result + 1.0) / 2.0;

  // Additional validation
  if (result < 0.0 || result > 1.0 || MathIsValidNumber (result) == false)
  {
    result = 0.2 + 0.6 * u.RNDprobab ();
  }

  return result;
}
//——————————————————————————————————————————————————————————————————————————————

TentMap 方法实现了帐篷映射，用于生成混沌序列。该方法接收输入值 x（要求在 0 到 1 之间），校验 x 的合法性；如有必要，用有效范围内的随机值替换。接下来，使用接近 2 的参数 mu，通过帐篷映射特有的分段线性函数计算出新的数值。

计算完成后，会再次校验数值有效性，必要时通过随机数进行归一化处理。随后方法返回新生成的混沌值。

//——————————————————————————————————————————————————————————————————————————————
double C_AO_COA_chaos::TentMap (double x)
{
  // Protection against incorrect inputs
  if (x < 0.0 || x > 1.0 || MathIsValidNumber (x) == false)
  {
    x = 0.2 + 0.6 * u.RNDprobab ();
  }

  // Tent map: x(n+1) = μ*min(x(n), 1-x(n))
  double mu = 1.99; // Parameter close to 2 for chaotic behavior
  double result;

  if (x <= 0.5) result = mu * x;
  else result = mu * (1.0 - x);

  // Additional validation
  if (result < 0.0 || result > 1.0 || MathIsValidNumber (result) == false)
  {
    result = 0.2 + 0.6 * u.RNDprobab ();
  }

  return result;
}
//——————————————————————————————————————————————————————————————————————————————

SelectChaosMap 方法用于根据指定类型选择并执行对应的混沌映射函数。该方法接收数值 x 和参数 type（用于指定具体的混沌映射类型）。方法的核心逻辑是：用 type 除以 3 的余数来选择映射方式，从而在三种混沌映射之间循环切换：逻辑映射、正弦映射和帐篷映射。根据计算结果，调用对应的函数，通过选定的混沌动力学规则将输入值 x 转换为新值。

如果因某种原因，type 超出预期范围（0、1、2），则默认使用逻辑映射。上述每一种映射都能模拟混沌行为，在优化过程中用于生成多样化、不可预测的随机数值。

//——————————————————————————————————————————————————————————————————————————————
double C_AO_COA_chaos::SelectChaosMap (double x, int type)
{
  // Select a chaotic map based on type
  switch (type % 3)
  {
    case 0:
      return LogisticMap (x);
    case 1:
      return SineMap (x);
    case 2:
      return TentMap (x);
    default:
      return LogisticMap (x);
  }
}
//——————————————————————————————————————————————————————————————————————————————

InitialPopulation 方法用于通过拉丁超立方抽样（LHS）技术初始化优化算法的初始种群。LHS是一种分层抽样方法，与纯随机抽样相比，它能更均匀地覆盖高维搜索空间，从而提升初始种群的质量。

方法首先声明两个数组：一个用于存储拉丁超立方数值，另一个用于辅助生成。方法会尝试为所需数组分配内存；若内存分配失败，则启用备用方案，通过随机方式生成初始种群。这一机制确保程序不会因内存不足而崩溃。

随后，方法生成拉丁超立方的抽样数值。针对每个坐标维度，先生成有序数值数组，再对其进行随机打乱。打乱后的数值被赋值给超立方数组。拉丁超立方的数值会被转换为搜索空间中个体的坐标。转换计算基于指定的区间范围完成，生成的数值会被限制在要求的取值范围和步长内。

最后，方法会设置一个标志位，标记初始种群已完成创建或修改。该方案的优势在于能够生成多样性更高的初始种群。

//——————————————————————————————————————————————————————————————————————————————
void C_AO_COA_chaos::InitialPopulation ()
{
  // Create Latin Hypercube for the initial population
  double latinCube []; // One-dimensional array for storing hypercube values
  double tempValues []; // Temporary array for storing and shuffling values

  ArrayResize (latinCube, popSize * coords);
  ArrayResize (tempValues, popSize);

  // Generate a Latin hypercube
  for (int c = 0; c < coords; c++)
  {
    // Create ordered values
    for (int i = 0; i < popSize; i++)
    {
      tempValues [i] = (double)i / popSize;
    }

    // Shuffle the values
    for (int i = popSize - 1; i > 0; i--)
    {
      int j = (int)(u.RNDprobab () * (i + 1));
      if (j < popSize)
      {
        double temp = tempValues [i];
        tempValues [i] = tempValues [j];
        tempValues [j] = temp;
      }
    }

    // Assign the mixed values
    for (int i = 0; i < popSize; i++)
    {
      latinCube [i * coords + c] = tempValues [i];
    }
  }

  // Convert the Latin hypercube values to coordinates
  for (int i = 0; i < popSize; i++)
  {
    for (int c = 0; c < coords; c++)
    {
      double x = rangeMin [c] + latinCube [i * coords + c] * (rangeMax [c] - rangeMin [c]);
      a [i].c [c] = u.SeInDiSp (x, rangeMin [c], rangeMax [c], rangeStep [c]);
    }
  }
}
//——————————————————————————————————————————————————————————————————————————————

FirstCarrierWaveSearch 方法实现了算法中的一个搜索阶段，旨在平衡全局空间探索与对已知优质解的局部开发。它的核心任务是更新智能体的位置与速度，持续挖掘并优化潜在解。方法开始时会定义一个系数，用于控制当前迭代轮次的探索强度。该系数会随着迭代轮次的增加呈二次递减，确保搜索重心从全局探索逐步转向局部优化。随后，针对种群中的每个智能体执行变异判断 —— 在设定概率下触发变异，以提升解的多样性。之后，对每一个搜索方向（坐标维度）执行以下操作：

• 选择用于生成新潜在解的混沌映射类型，
• 选择全局搜索或局部搜索策略。

在全局搜索模式下，智能体通过混沌分量更新位置，并在调整速度时兼顾惯性与运动方向。在局部搜索模式下，智能体聚焦于已找到的最优解，向其进行加权靠拢，并加入少量随机扰动以避免循环停滞。两种模式下都会对速度进行限制，防止智能体大幅跳出搜索空间边界。智能体的位置会根据搜索空间约束进行更新；若检测到约束违反，则执行修正操作。此时会调整位置并降低速度，使后续的移动步长更加平缓。

//——————————————————————————————————————————————————————————————————————————————
void C_AO_COA_chaos::FirstCarrierWaveSearch ()
{
  // Adaptive balance between exploration and exploitation
  double globalPhase = (double)epochNow / S1;
  double explorationRate = 1.0 - globalPhase * globalPhase; // Quadratic decrease

  // For each agent
  for (int i = 0; i < popSize; i++)
  {
    // Apply mutations with some probability to increase diversity
    if (u.RNDprobab () < mutationRate * (1.0 + explorationRate))
    {
      ApplyMutation (i);
      continue;
    }

    for (int c = 0; c < coords; c++)
    {
      // Select a chaotic map with uniform distribution
      int mapType = ((i + c + epochNow) % 3);

      // Safely check access to the gamma array
      if (c < ArraySize (agent [i].gamma))
      {
        agent [i].gamma [c] = SelectChaosMap (agent [i].gamma [c], mapType);
      }
      else
      {
        continue; // Skip if the index is invalid
      }

      // Determine the relationship between global and local search
      double strategy = u.RNDprobab ();
      double x;

      if (strategy < explorationRate)
      {
        // Global search with a chaotic component
        x = rangeMin [c] + agent [i].gamma [c] * (rangeMax [c] - rangeMin [c]);

        // Add a velocity component to maintain the movement direction
        agent [i].velocity [c] = inertia * agent [i].velocity [c] +
                                 (1.0 - inertia) * (x - a [i].c [c]);
      }
      else
      {
        // Local search around the best solutions
        double personalAttraction = u.RNDprobab ();
        double globalAttraction = u.RNDprobab ();

        // Balanced attraction to the best solutions
        double attractionTerm = //personalAttraction * (agent [i].cPrev [c] - a [i].c [c]) +
                                personalAttraction * (a [i].cB [c] - a [i].c [c]) +
                                socialFactor * globalAttraction * (cB [c] - a [i].c [c]);

        // Chaotic disturbance to prevent being stuck
        double chaosRange = alpha [c] * explorationRate;
        double chaosTerm = chaosRange * (2.0 * agent [i].gamma [c] - 1.0);

        // Update velocity with inertia
        agent [i].velocity [c] = inertia * agent [i].velocity [c] +
                                 (1.0 - inertia) * (attractionTerm + chaosTerm);
      }

      // Limit the velocity to prevent too large steps
      double maxVelocity = 0.1 * (rangeMax [c] - rangeMin [c]);
      if (MathAbs (agent [i].velocity [c]) > maxVelocity)
      {
        agent [i].velocity [c] = maxVelocity * (agent [i].velocity [c] > 0 ? 1.0 : -1.0);
      }

      // Apply the velocity to the position
      x = a [i].c [c] + agent [i].velocity [c];

      // Apply search space restrictions
      a [i].c [c] = u.SeInDiSp (x, rangeMin [c], rangeMax [c], rangeStep [c]);

      // Check the constraints and apply a smooth correction
      double violation = CalculateConstraintValue (i, c);
      if (violation > eps)
      {
        double gradient = CalculateWeightedGradient (i, c);
        double correction = -gradient * violation * (1.0 - globalPhase);
        a [i].c [c] = u.SeInDiSp (a [i].c [c] + correction, rangeMin [c], rangeMax [c], rangeStep [c]);

        // Reset the velocity when correcting violations
        agent [i].velocity [c] *= 0.5;
      }
    }
  }
}
//——————————————————————————————————————————————————————————————————————————————

SecondCarrierWaveSearch 方法是初始搜索之后的优化阶段，用于对已找到的解进行深化与精细化寻优。该方法的核心目标，是通过更精细的搜索策略与参数自适应，提升上一阶段得到的优化结果。

方法首先计算一个反映局部搜索阶段的参数，该参数会随迭代逐步增强。这使得算法从大范围搜索，逐步转向对已知优质解所在区域的更细致、更精准的探索。方法开头会先检查算法是否已收敛。如果算法进入稳定状态，会对部分智能体执行变异，以提升解的多样性，避免陷入局部最优。

对每个智能体，在其搜索空间内进行逐维度的新解搜寻。算法会选定混沌映射，用于引入随机扰动。搜索参数会随着逼近最优解而逐步减小。从而在当前最优解附近实现更聚焦的局部搜索。在更新每一个位置时，算法会保留智能体此前的寻优成果。系统会确定一个基准点，该基准点可以是全局最优解，也可以是智能体自身的历史最优解，从而兼顾个体经验与种群整体成果。

位置更新过程同时使用混沌偏移与随机噪声（如莱维噪声），进一步增强随机性，有助于发现更优的新解。方法在更新智能体速度时会引入惯性，使位置变化更平滑，避免剧烈跳跃。最终，更新后的位置会被约束在设定范围内，保证满足问题的约束条件。

简言之，SecondCarrierWaveSearch 方法致力于对已有解进行更精准、更深入的优化。

//——————————————————————————————————————————————————————————————————————————————
void C_AO_COA_chaos::SecondCarrierWaveSearch ()
{
  // Refining local search with adaptive parameters
  double localPhase = (double)(epochNow - S1) / S2;
  double intensificationRate = localPhase * localPhase; // Quadratic increase in intensification

  // Check the algorithm convergence
  bool isConverged = IsConverged ();

  // For each agent
  for (int i = 0; i < popSize; i++)
  {
    // If convergence is detected, add a random mutation to some agents
    if (isConverged && i % 3 == 0)
    {
      ApplyMutation (i);
      continue;
    }

    for (int c = 0; c < coords; c++)
    {
      // Select a chaotic map with uniform distribution
      int mapType = ((i * c + epochNow) % 3);
      agent [i].gamma [c] = SelectChaosMap (agent [i].gamma [c], mapType);

      // Adaptive search radius with narrowing towards the end of optimization
      double adaptiveAlpha = alpha [c] * (1.0 - 0.8 * intensificationRate);

      // Select a base point with priority to the best solutions
      double basePoint;
      if (a [i].f > a [i].fB)
      {
        basePoint = a [i].c [c];  // The current position is better
      }
      else
      {
        double r = u.RNDprobab ();

        if (r < 0.7 * (1.0 + intensificationRate)) // Increase attraction to the global best
        {
          basePoint = cB [c];  // Global best
        }
        else
        {
          basePoint = a [i].cB [c];  // Personal best
        }
      }

      // Local search with a chaotic component
      double chaosOffset = adaptiveAlpha * (2.0 * agent [i].gamma [c] - 1.0);

      // Add Levy noise for random long jumps (heavy tailed distribution) 
      double levyNoise = 0.0;
      if (u.RNDprobab () < 0.1 * (1.0 - intensificationRate))
      {
        // Simplified approximation of Levy noise
        double u1 = u.RNDprobab ();
        double u2 = u.RNDprobab ();

        if (u2 > 0.01) // Protection against division by very small numbers
        {
          levyNoise = 0.01 * u1 / MathPow (u2, 0.5) * adaptiveAlpha * (rangeMax [c] - rangeMin [c]);
        }
      }

      // Update the velocity with inertia
      agent [i].velocity [c] = inertia * (1.0 - 0.5 * intensificationRate) * agent [i].velocity [c] +
                               (1.0 - inertia) * (chaosOffset + levyNoise);

      // Apply the velocity to the position
      double x = basePoint + agent [i].velocity [c];

      // Limit the position
      a [i].c [c] = u.SeInDiSp (x, rangeMin [c], rangeMax [c], rangeStep [c]);
    }
  }
}
//——————————————————————————————————————————————————————————————————————————————

在下一篇文章中，我们将继续讲解该算法的其余方法，进行测试并总结结果。