2022高教社杯全国大学生数学建模竞赛A题第二问第一小问解析+代码

一、引言

2022年A题是围绕波浪能发电装置进行运动学建模，整体建模难度中等；可以调用模型进行求解，编程难度中等，整体解题难度中等偏下。

第二问紧接着第一问，根据第一问的建模进行优化，可以调用多种优化模型进行优化。第二问分为两个部分，阻尼力不变和阻尼力变化。阻尼力不变时候为一维线性优化，阻尼力变化为双变量二维优化。

题目下载链接：2022高教社杯全国大学生数学建模竞赛赛题下载(官方链接)

二、建模思路

第一步：物理系统与假设

• 研究对象：与问题1相同，为一个波浪能转换装置，由外部浮子、内部振子以及连接二者的能量输出系统（PTO）构成。PTO包含弹簧和直线阻尼器。
• 核心假设：
  1. 浮子只做垂荡运动，就是纯上下平移。
  2. 海水无粘、无旋，波浪为线性微幅波。
  3. 忽略所有摩擦及PTO等部件的质量。
  4. 核心变化：波浪频率更换为 $\omega =2.2143{\text{s}}^{-1}$，目标是优化阻尼系数以最大化平均输出功率。

第二步：变量与参数定义（同第一问）

• 状态变量：
  • 浮子垂荡位移和速度：$x_f(t)$, $v_f(t)$
  • 振子垂荡位移和速度：$x_o(t)$, $v_o(t)$
• 关键参数：
  • 系统参数：$m_f,m_o,m_a,k_s,c_w,B,F_0$。
  • 优化变量：直线阻尼器的阻尼系数 $c_d$。本问题旨在寻找使平均功率最大的 $c_d$ 值。

第三步：受力分析与功率定义

浮子与振子的受力分析与问题1完全相同。因此，我们直接得到运动方程。

功率计算是问题2的新核心：
阻尼器是能量提取元件，其瞬时功率等于阻尼力乘以相对速度：
$$P(t)=F_d(t)\cdot v_{\text{rel}}(t)$$
其中，$v_{\text{rel}}(t)=v_o(t)-v_f(t)$。

• 对于情形1（线性阻尼）：$F_d=c_d\cdot v_{\text{rel}}$，代入得：
  $$P(t)=c_d\cdot (v_{\text{rel}}(t){)}^2$$
• 优化目标是最大化平均输出功率，即在一段稳定时间 $T$ 内对瞬时功率取平均：
  $$\bar{P}=\frac{1}{T}{\int }_{t_0}^{t_0+T}P(t)dt\approx \frac{1}{N}\sum _{i=1}^N{c}_d\cdot (v_{\text{rel}}(t_i){)}^2$$

第四步：运动方程与优化模型建立

1. 运动微分方程（线性阻尼）：

   浮子方程：
   $$(m_f+m_a){\ddot{x}}_f=F_0\cos (\omega t)+c_d(v_o-v_f)+k_s(x_o-x_f)-c_w{v}_f-B{x}_f$$

   振子方程：
   $$m_o{\ddot{x}}_o=-c_d(v_o-v_f)-k_s(x_o-x_f)$$

2. 优化模型建立：

   • 优化变量：$c_d$
   • 目标函数：平均输出功率 $\bar{P}(c_d)$
   • 约束条件：$0\le c_d\le 100000$，且运动满足上述微分方程及初始条件。
   • 问题表述：寻找一个最优的 $c_d^{\ast }$，使得系统在给定波浪激励下，从运动方程解出的 $\bar{P}(c_d)$ 达到最大。
     $$c_d^{\ast }=\arg \underset{0\le c_d\le 10^5}{\max }\bar{P}(c_d)$$

第五步：求解策略

1. 为什么需要数值优化？
   目标函数 $\bar{P}(c_d)$ 是隐式的。你无法直接写出一个关于 $c_d$ 的简单函数式，因为 $\bar{P}$ 的值依赖于微分方程组的数值解。这构成了一个仿真优化问题。

2. 求解的基本思想（两层循环）：

   • 内层循环（仿真）：对于一个给定的候选阻尼系数 $c_d$，通过龙格-库塔法求解微分方程组，得到系统在足够长时间（如100s后达到稳态）内的运动响应 $x_f(t),x_o(t),v_f(t),v_o(t)$。
   • 外层循环（优化）：采用优化算法（如遗传算法，如你所提供的代码中使用的）在区间 $[0,100000]$ 内搜索 $c_d$。该算法会生成一系列 $c_d$ 的候选值，对每个值调用内层循环进行仿真并计算其对应的平均功率 $\bar{P}$，然后根据功率大小引导搜索方向，最终逼近使 $\bar{P}$ 最大的 $c_d^{\ast }$。

三、编程计算

完整代码下载

遗传算法代码

    def objective_function(self, individual: List[float]) -> Tuple[float]:
        """目标函数 - 计算平均功率（最大化问题）"""
        try:
            power = self.run_simulation(individual[0])
            return (power,)
        except Exception as e:
            print(f"模拟错误: {e}")
            return (0.0,)
    
    def setup_genetic_algorithm(self) -> base.Toolbox:
        """设置遗传算法工具箱"""
        creator.create("FitnessMax", base.Fitness, weights=(1.0,))
        creator.create("Individual", list, fitness=creator.FitnessMax)
        
        toolbox = base.Toolbox()
        toolbox.register("attr_float", random.uniform, 0, 100000)
        toolbox.register("individual", tools.initRepeat, creator.Individual, 
                         toolbox.attr_float, n=1)
        toolbox.register("population", tools.initRepeat, list, toolbox.individual)
        toolbox.register("evaluate", self.objective_function)
        toolbox.register("mate", tools.cxBlend, alpha=0.5)
        toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1000, indpb=0.2)
        toolbox.register("select", tools.selTournament, tournsize=3)
        
        return toolbox