1. 项目概述当医学影像遇上能源优化看到“2026年河北省研究生数学建模 C/D 题”这个标题很多同学的第一反应可能是“又来画饼了”。确实预测两年后的赛题听起来有点玄乎但作为一名带过好几届数模竞赛、自己也从坑里爬出来的老手我得说这个标题背后透露出的命题趋势和备赛逻辑恰恰是现在准备竞赛最需要关注的核心。它把看似风马牛不相及的“冠脉造影狭窄识别”和“零碳园区微电网优化”放在一起绝不是随意拼凑而是精准地指向了当前交叉学科应用和复杂系统建模的两大前沿热点。简单来说这道“预测题”模拟了研究生数模竞赛的典型风格一道偏向数据科学与人工智能的C题一道偏向运筹优化与系统控制的D题。C题要求你从医学影像中自动识别血管狭窄这本质是一个计算机视觉中的图像分割与分类问题考验的是数据预处理、特征工程和模型构建能力。而D题则要求你对一个包含光伏、风电、储能和负荷的园区微电网进行优化调度目标是实现零碳运行与经济性的平衡这涉及到混合整数规划、动态优化甚至强化学习等运筹学方法。两者共同的核心都是将现实世界的复杂问题抽象为可计算的数学模型并给出可验证的解决方案。无论你是计算机、电气、自动化专业还是数学、统计背景这道题都能找到你的发力点。对于准备参赛的同学提前琢磨这类题目不是为了押题而是为了掌握一种“解题框架思维”——即面对一个陌生领域的具体问题如何快速进行问题分解、方法选型和模型构建。接下来我就结合自己多年踩坑和指导的经验对这两个题目进行深度拆解并提供一套可直接上手、能灵活调整的代码框架与建模思路。2. 核心思路拆解从问题到模型的映射逻辑面对一个综合性的数模题目最忌讳的就是一头扎进细节里开始编程。正确的打开方式是先进行顶层设计明确从“问题描述”到“最终论文”每一步要做什么用什么方法做以及为什么这么做。2.1 C题冠脉造影狭窄识别——一个标准的CV Pipeline医学影像分析题目在近年的国赛、美赛中屡见不鲜。冠脉造影狭窄识别的核心任务可以拆解为输入一张冠脉造影的X光影像输出图像中血管的轮廓以及狭窄部位的位置和程度评估。这本质上是一个“语义分割”“关键点/区域分类”的组合任务。为什么选择语义分割而不是单纯的目标检测因为血管是细长的管状结构狭窄表现为血管局部内径的变细。仅仅用检测框标出狭窄区域目标检测无法量化狭窄程度比如直径狭窄率。我们需要得到血管像素级的精确轮廓语义分割然后基于轮廓计算血管各处的直径从而找出直径小于正常值一定比例通常临床以50%或70%为界的狭窄段。整体建模Pipeline设计如下数据预处理与增强医学影像质量参差不齐存在对比度低、噪声大、造影剂不均匀等问题。预处理是关键第一步。血管分割模型采用编码器-解码器结构的神经网络如U-Net及其变体进行像素级分割将图像分为“血管”和“背景”两类。血管中心线提取与直径分析对分割出的二值血管掩膜图像采用骨架化算法提取单像素宽的中心线。沿着中心线法线方向测量血管的局部直径。狭窄识别与分类根据血管直径剖面识别局部最小值。计算狭窄率(1 - 狭窄处直径 / 参考正常直径) * 100%。根据狭窄率阈值如50%为轻度50%-70%为中度70%为重度对狭窄部位进行分类。结果可视化与评估将分割的血管轮廓、狭窄定位框和分类结果以热力图或标注框形式叠加在原图上并计算分割的Dice系数、狭窄识别的准确率/召回率等指标。注意临床实际中参考正常直径的选取有讲究。通常选取狭窄近端和远端看起来正常的血管段直径的平均值或采用同一血管的未狭窄段直径。在建模时如果数据未标注正常段可以假设血管直径变化平滑采用滑动窗口内的局部最大直径或拟合的基线直径作为参考。2.2 D题零碳园区微电网优化——一个典型的混合整数规划问题零碳园区微电网优化是一个多时间尺度、多设备耦合、多目标约束的优化问题。核心目标是在满足园区电负荷需求的前提下协调调度光伏发电、风力发电、蓄电池储能以及可能的外部电网购电假设外部电网为高碳电网购电会增加碳足迹使得一个调度周期如24小时内的总运行成本最低同时实现净零碳排放。为什么这是一个混合整数规划MIP问题因为模型中既包含连续变量如发电功率、储能充放电功率、荷电状态也包含整数变量或0-1变量。典型的整数变量场景包括储能的充放电状态同一时刻不能既充又放需要引入0-1变量表示模式、燃气轮机等可启停机组的启停状态、与电网的买卖电状态防止同时买卖等。这些“要么开要么关”的决策需要用整数变量来刻画。模型核心要素拆解决策变量各时刻光伏/风电出力通常视为不可控但可预测、储能充/放电功率、储能荷电状态SOC、从电网购电功率、向电网售电功率如果允许、可能的备用容量等。目标函数最小化总成本 购电成本 - 售电收入 设备运维成本与充放电量相关 储能折旧成本。为了实现“零碳”通常将碳排放量作为约束条件而非目标函数。约束条件功率平衡约束每一时刻发电总量 放电功率 购电功率 负荷 充电功率 售电功率。设备运行约束光伏/风电出力上下限储能SOC上下限、充放电功率上下限、充放电效率、初始/终止SOC约束充放电互斥约束需0-1变量。电网交互约束购电/售电功率上下限可能存在的功率因数约束。零碳约束调度周期内总碳排放量 ≤ 0。碳排放来源主要是外购电力假设其有碳强度系数如kgCO₂/kWh。园区内光伏、风电、储能视为零碳。因此约束可转化为外购电量 * 碳强度系数 ≤ 0这实质上强制要求外购电量为零除非题目允许碳交易即购买绿证或碳汇来抵消。这是一个非常强的约束会极大改变优化策略。不确定性处理光伏、风电和负荷预测都存在误差。高级的模型会考虑鲁棒优化或随机规划引入不确定性集合或场景使调度方案对预测误差不敏感。对于初次接触可以先做确定性优化。3. C题代码框架与实操要点这里提供一个基于Python和PyTorch的简化版代码框架重点展示核心流程和易错点。3.1 数据准备与预处理模块医学影像数据通常以DICOM格式存储但竞赛为了方便很可能提供的是处理后的JPEG或PNG图像甚至可能是公开数据集如部分标注的冠脉造影数据集。import cv2 import numpy as np import torch from torch.utils.data import Dataset, DataLoader import albumentations as A from albumentations.pytorch import ToTensorV2 class CoronaryAngioDataset(Dataset): def __init__(self, image_paths, mask_pathsNone, is_trainTrue, img_size512): self.image_paths image_paths self.mask_paths mask_paths self.is_train is_train self.img_size img_size # 训练和验证/测试采用不同的数据增强策略 if self.is_train: self.transform A.Compose([ A.Resize(heightimg_size, widthimg_size), A.RandomRotate90(p0.5), A.HorizontalFlip(p0.5), A.VerticalFlip(p0.5), # 医学影像增强对比度受限自适应直方图均衡化(CLAHE)非常有效 A.CLAHE(clip_limit2.0, tile_grid_size(8,8), p0.8), A.RandomBrightnessContrast(brightness_limit0.1, contrast_limit0.1, p0.5), A.GaussNoise(var_limit(10.0, 50.0), p0.3), A.Normalize(mean[0.485], std[0.229]), # 灰度图单通道 ToTensorV2(), ]) else: self.transform A.Compose([ A.Resize(heightimg_size, widthimg_size), A.CLAHE(clip_limit2.0, tile_grid_size(8,8), always_applyTrue), # 测试时也做CLAHE A.Normalize(mean[0.485], std[0.229]), ToTensorV2(), ]) def __getitem__(self, idx): image cv2.imread(self.image_paths[idx], cv2.IMREAD_GRAYSCALE) if image is None: raise FileNotFoundError(f无法读取图像: {self.image_paths[idx]}) # 确保图像是二维的 if len(image.shape) 3: image cv2.cvtColor(image, cv2.COLOR_BGR2GRAY) if self.mask_paths is not None: mask cv2.imread(self.mask_paths[idx], cv2.IMREAD_GRAYSCALE) mask (mask 127).astype(np.uint8) # 二值化 augmented self.transform(imageimage, maskmask) image_tensor augmented[image] mask_tensor augmented[mask].float().unsqueeze(0) # 增加通道维 return image_tensor, mask_tensor else: augmented self.transform(imageimage) return augmented[image], self.image_paths[idx] # 测试时返回路径用于保存结果 def __len__(self): return len(self.image_paths)实操心得CLAHE是关键对于对比度低的X光影像CLAHE比简单的直方图均衡化效果更好能增强血管与背景的对比且不易过度增强噪声。谨慎使用几何变换旋转、翻转对于血管结构是有效的因为血管方向不固定。但像弹性形变这类强变换要小心可能扭曲血管的拓扑结构如分叉点。二值化阈值标注的掩膜可能是0-255的灰度图需要统一二值化。127是个常用阈值但最好检查一下标注数据的实际范围。3.2 模型构建U-Net及其改进U-Net是医学分割的经典网络。这里实现一个基础版并讨论几个改进方向。import torch.nn as nn import torch.nn.functional as F class DoubleConv(nn.Module): (卷积 BN ReLU) * 2 def __init__(self, in_channels, out_channels): super().__init__() self.double_conv nn.Sequential( nn.Conv2d(in_channels, out_channels, kernel_size3, padding1), nn.BatchNorm2d(out_channels), nn.ReLU(inplaceTrue), nn.Conv2d(out_channels, out_channels, kernel_size3, padding1), nn.BatchNorm2d(out_channels), nn.ReLU(inplaceTrue) ) def forward(self, x): return self.double_conv(x) class UNet(nn.Module): def __init__(self, n_channels1, n_classes1): super(UNet, self).__init__() self.n_channels n_channels self.n_classes n_classes self.inc DoubleConv(n_channels, 64) self.down1 nn.Sequential(nn.MaxPool2d(2), DoubleConv(64, 128)) self.down2 nn.Sequential(nn.MaxPool2d(2), DoubleConv(128, 256)) self.down3 nn.Sequential(nn.MaxPool2d(2), DoubleConv(256, 512)) self.down4 nn.Sequential(nn.MaxPool2d(2), DoubleConv(512, 1024)) self.up1 nn.ConvTranspose2d(1024, 512, kernel_size2, stride2) self.conv1 DoubleConv(1024, 512) # 1024 512(skip) 512(up) self.up2 nn.ConvTranspose2d(512, 256, kernel_size2, stride2) self.conv2 DoubleConv(512, 256) self.up3 nn.ConvTranspose2d(256, 128, kernel_size2, stride2) self.conv3 DoubleConv(256, 128) self.up4 nn.ConvTranspose2d(128, 64, kernel_size2, stride2) self.conv4 DoubleConv(128, 64) self.outc nn.Conv2d(64, n_classes, kernel_size1) def forward(self, x): x1 self.inc(x) x2 self.down1(x1) x3 self.down2(x2) x4 self.down3(x3) x5 self.down4(x4) x self.up1(x5) # 拼接跳跃连接注意空间尺寸要对齐 x torch.cat([x, x4], dim1) x self.conv1(x) x self.up2(x) x torch.cat([x, x3], dim1) x self.conv2(x) x self.up3(x) x torch.cat([x, x2], dim1) x self.conv3(x) x self.up4(x) x torch.cat([x, x1], dim1) x self.conv4(x) logits self.outc(x) return logits # 改进方向讨论 # 1. 深度监督在中间层也添加辅助输出帮助梯度回传缓解梯度消失。 # 2. 注意力机制在跳跃连接处加入注意力门Attention Gate让网络更关注血管区域抑制背景噪声。 # 3. 使用预训练编码器将U-Net的编码器替换为ResNet、EfficientNet等在大规模图像数据上预训练的网络可以提升特征提取能力尤其在小数据集上效果显著。 # 4. 损失函数二值分割常用Dice Loss BCE Loss的组合。Dice Loss直接优化分割区域的重叠度对类别不平衡血管像素少问题更鲁棒。3.3 后处理与狭窄分析模块这是将分割结果转化为最终答案的关键也是容易出bug的地方。import cv2 import numpy as np from skimage import morphology, measure def analyze_stenosis_from_mask(binary_mask): 从二值血管掩膜分析狭窄情况。 参数: binary_mask: numpy数组二值图像1为血管0为背景。 返回: results: 列表每个元素是一个字典描述一个狭窄区域。 # 1. 细化骨架化获取中心线 # 使用中轴骨架化保留拓扑结构 skeleton morphology.skeletonize(binary_mask 0) # 2. 获取中心线的坐标点 skeleton_points np.column_stack(np.where(skeleton 0)) # 3. 构建中心线的图结构用于追踪血管分支 # 这里简化处理假设血管是单连通的无分叉或处理主干。 # 对于复杂分叉需要先对骨架进行端点、交叉点分析分割成多个段。 from scipy.spatial import KDTree tree KDTree(skeleton_points) # 4. 沿着中心线测量局部直径 diameters [] for i, point in enumerate(skeleton_points): # 对于中心线上每个点计算到最近背景点的距离*2 作为该点处血管直径的近似 # 更准确的做法在该点法线方向需估算骨架方向上向两边搜索背景点。 # 这里使用距离变换作为简化方法。 pass # 具体实现需要结合距离变换 # 5. 平滑直径曲线并寻找局部极小值狭窄点 from scipy.signal import savgol_filter, find_peaks smoothed_diameters savgol_filter(diameters, window_length11, polyorder3) # 寻找局部最小值狭窄注意find_peaks找的是极大值对负值取反找极小值 minima_indices, _ find_peaks(-smoothed_diameters, prominence0.5) # prominence需要根据数据调整 # 6. 计算狭窄率并分类 results [] for idx in minima_indices: stenosis_diameter smoothed_diameters[idx] # 参考直径取狭窄点前后一段正常血管直径的平均值例如前后各10个点排除其他狭窄点 window_size 10 start max(0, idx - window_size) end min(len(smoothed_diameters), idx window_size 1) # 简单起见取窗口内最大值作为参考直径假设血管直径变化平缓 reference_diameter np.max(smoothed_diameters[start:end]) stenosis_ratio (1 - stenosis_diameter / reference_diameter) * 100 severity 轻度 if stenosis_ratio 70: severity 重度 elif stenosis_ratio 50: severity 中度 # 获取狭窄点在原图中的坐标近似 stenosis_point skeleton_points[idx] results.append({ location_xy: (int(stenosis_point[1]), int(stenosis_point[0])), # (x, y) stenosis_ratio_percent: round(stenosis_ratio, 2), severity: severity, reference_diameter: round(reference_diameter, 2), stenosis_diameter: round(stenosis_diameter, 2) }) return results # 可视化函数 def visualize_results(original_image, binary_mask, stenosis_results, save_pathNone): vis_img cv2.cvtColor(original_image, cv2.COLOR_GRAY2BGR) # 用半透明红色覆盖血管区域 mask_color np.zeros_like(vis_img) mask_color[binary_mask 0] [0, 0, 255] # BGR中的红色 cv2.addWeighted(mask_color, 0.3, vis_img, 0.7, 0, vis_img) # 标记狭窄点 for res in stenosis_results: x, y res[location_xy] cv2.circle(vis_img, (x, y), radius5, color(0, 255, 0), thickness-1) # 绿色圆点 label f{res[severity]}({res[stenosis_ratio_percent]}%) cv2.putText(vis_img, label, (x10, y-10), cv2.FONT_HERSHEY_SIMPLEX, 0.5, (0, 255, 0), 1) if save_path: cv2.imwrite(save_path, vis_img) return vis_img踩坑实录骨架化噪声二值掩膜边缘如果有毛刺骨架化会产生大量短小分支干扰中心线提取。务必先对掩膜进行morphology.opening开运算去除小噪声。血管分叉处理真实的冠脉血管是有分叉的。简单的骨架化会在分叉点产生一个“星形”结构导致中心线追踪混乱。需要在骨架化后进行端点、交叉点检测将血管分解成多个独立的“段”再分别分析。可以使用skimage.morphology中的remove_small_objects和skan库Skeleton Analysis进行更专业的骨架分析。直径测量不准在弯曲的血管处法线方向估算不准会导致直径测量偏差。一个更稳健的方法是使用距离变换dist_transform cv2.distanceTransform(binary_mask.astype(np.uint8), cv2.DIST_L2, 3)。对于中心线上的点其距离变换值的两倍可以较好地近似该点处的血管半径。但要注意在血管边缘不平滑时仍有误差。狭窄判定阈值临床上的50%、70%阈值是金标准但在模型评估时如果数据标注的“狭窄”是基于视觉评估的可能与计算的直径狭窄率有系统偏差。建议在验证集上做一个简单的校准看看模型计算的狭窄率与标注的狭窄等级如果有的话之间的对应关系微调阈值。4. D题建模与求解框架微电网优化建模通常使用线性规划LP或混合整数线性规划MILP工具。Python中PuLP和CVXPY对于线性问题是易上手的建模库而Gurobi、CPLEX是强大的商业求解器学术可免费申请license。这里以PuLP为例构建一个确定性优化模型。4.1 模型参数与变量定义假设调度周期T24小时时间间隔Δt1小时。import pulp import numpy as np # 参数定义 (示例值需根据题目数据修改) T 24 # 时间步数 # 预测数据 (kW) P_pv_pred np.array([...]) # 光伏预测出力长度T P_wind_pred np.array([...]) # 风电预测出力 P_load_pred np.array([...]) # 电负荷预测 # 电价 (元/kWh) price_buy np.array([...]) # 购电价分时电价 price_sell np.array([...]) # 售电价通常低于购电价 carbon_intensity_grid 0.8 # 电网碳排放强度 kgCO2/kWh # 储能参数 E_batt_max 500 # 储能最大容量 kWh E_batt_min 50 # 储能最小容量防止过放设为最大容量的10% P_charge_max 100 # 最大充电功率 kW P_discharge_max 100 # 最大放电功率 kW eta_charge 0.95 # 充电效率 eta_discharge 0.95 # 放电效率 E_batt_initial 250 # 初始电量 kWh E_batt_end 250 # 调度周期末期望电量 kWh通常等于初始值 # 创建问题 prob pulp.LpProblem(ZeroCarbon_Microgrid_Optimization, pulp.LpMinimize) # 决策变量 # 连续变量 P_grid_buy pulp.LpVariable.dicts(P_grid_buy, range(T), lowBound0) # 购电功率 P_grid_sell pulp.LpVariable.dicts(P_grid_sell, range(T), lowBound0) # 售电功率 P_batt_charge pulp.LpVariable.dicts(P_batt_charge, range(T), lowBound0, upBoundP_charge_max) P_batt_discharge pulp.LpVariable.dicts(P_batt_discharge, range(T), lowBound0, upBoundP_discharge_max) E_batt pulp.LpVariable.dicts(E_batt, range(T1), lowBoundE_batt_min, upBoundE_batt_max) # 从t0到tT # 0-1整数变量 (确保不同时充放电) u_charge pulp.LpVariable.dicts(u_charge, range(T), catBinary) # 1表示充电 u_discharge pulp.LpVariable.dicts(u_discharge, range(T), catBinary) # 1表示放电 # 目标函数最小化总运行成本 # 成本 购电成本 - 售电收入 储能循环损耗成本简化为与充放电量成正比 cost_buy pulp.lpSum([P_grid_buy[t] * price_buy[t] for t in range(T)]) cost_sell pulp.lpSum([P_grid_sell[t] * price_sell[t] for t in range(T)]) # 假设储能充放电损耗成本系数为0.05元/kWh cost_batt_cycle 0.05 * pulp.lpSum([P_batt_charge[t] P_batt_discharge[t] for t in range(T)]) prob cost_buy - cost_sell cost_batt_cycle4.2 约束条件构建这是模型的核心需要严谨地表达所有物理和运行限制。# 约束条件 for t in range(T): # 1. 功率平衡约束 (核心) prob (P_pv_pred[t] P_wind_pred[t] P_batt_discharge[t] P_grid_buy[t] P_load_pred[t] P_batt_charge[t] P_grid_sell[t]), fPower_Balance_t{t} # 2. 储能电量动态方程 if t 0: prob (E_batt[t] E_batt_initial), fBatt_Init prob (E_batt[t1] E_batt[t] P_batt_charge[t] * eta_charge * 1 # Δt1小时能量功率*时间 - P_batt_discharge[t] / eta_discharge * 1), fBatt_Dynamics_t{t} # 3. 充放电互斥约束 (Big-M法) M max(P_charge_max, P_discharge_max) * 1.5 # 一个足够大的数 prob (P_batt_charge[t] u_charge[t] * M), fCharge_Logic_1_t{t} prob (P_batt_discharge[t] u_discharge[t] * M), fDischarge_Logic_1_t{t} prob (u_charge[t] u_discharge[t] 1), fCharge_Discharge_Mutex_t{t} # 4. 电网交互约束 (假设最大交互功率为200kW) P_grid_max 200 prob (P_grid_buy[t] P_grid_max), fGrid_Buy_Max_t{t} prob (P_grid_sell[t] P_grid_max), fGrid_Sell_Max_t{t} # 通常不允许同时购售电但上述互斥约束已通过功率平衡和成本函数自然实现因为售电价格低同时购售电不经济。如需严格约束可再加0-1变量。 # 储能周期首尾电量约束保证可持续运行 prob (E_batt[T] E_batt_end), Batt_End_Min # 也可设为等于强制要求周期末回到初始状态 # prob (E_batt[T] E_batt_end), Batt_End_Equal # 零碳约束 # 总碳排放 总购电量 * 碳强度 total_carbon_emission pulp.lpSum([P_grid_buy[t] for t in range(T)]) * carbon_intensity_grid # 强制零碳总碳排放 0 prob (total_carbon_emission 0), Zero_Carbon_Constraint # 注意这个约束非常强意味着 P_grid_buy[t] 对所有t都必须为0。 # 如果题目允许少量碳排放或碳交易约束应改为 Carbon_Budget (碳排放预算)。4.3 模型求解与结果分析# 求解 # 指定求解器例如CBC开源或Gurobi需安装且配置license solver pulp.PULP_CBC_CMD(msgFalse) # 使用CBC不输出求解日志 # solver pulp.GUROBI_CMD(msgTrue) # 如果安装了Gurobi prob.solve(solver) # 输出结果 print(f求解状态: {pulp.LpStatus[prob.status]}) print(f最优总成本: {pulp.value(prob.objective):.2f} 元) if prob.status pulp.LpOptimal: # 提取变量值 P_buy_opt [pulp.value(P_grid_buy[t]) for t in range(T)] P_sell_opt [pulp.value(P_grid_sell[t]) for t in range(T)] P_charge_opt [pulp.value(P_batt_charge[t]) for t in range(T)] P_discharge_opt [pulp.value(P_batt_discharge[t]) for t in range(T)] E_batt_opt [pulp.value(E_batt[t]) for t in range(T1)] # 计算实际碳排放 total_buy sum(P_buy_opt) actual_carbon total_buy * carbon_intensity_grid print(f总购电量: {total_buy:.2f} kWh, 理论碳排放: {actual_carbon:.2f} kgCO2) # 在零碳约束下total_buy应为0 # 可视化调度计划 import matplotlib.pyplot as plt fig, axes plt.subplots(3, 1, figsize(12, 10)) t_range list(range(T)) axes[0].plot(t_range, P_load_pred, label负荷, colorblack, linewidth2) axes[0].plot(t_range, P_pv_pred, label光伏, colorgold) axes[0].plot(t_range, P_wind_pred, label风电, colorcyan) axes[0].set_ylabel(功率 (kW)) axes[0].set_title(预测的发电与负荷) axes[0].legend() axes[0].grid(True, alpha0.3) axes[1].bar(t_range, P_charge_opt, label充电, colorgreen, alpha0.6) axes[1].bar(t_range, [-x for x in P_discharge_opt], label放电, colorred, alpha0.6) # 放电为负 axes[1].plot(t_range, E_batt_opt[:-1], label储能电量, colorblue, markero) axes[1].set_ylabel(功率/电量 (kW/kWh)) axes[1].set_title(储能调度) axes[1].legend() axes[1].grid(True, alpha0.3) axes[2].bar(t_range, P_buy_opt, label购电, colororange, alpha0.6) axes[2].bar(t_range, [-x for x in P_sell_opt], label售电, colorpurple, alpha0.6) # 售电为负 axes[2].set_xlabel(时间 (h)) axes[2].set_ylabel(功率 (kW)) axes[2].set_title(与电网交互) axes[2].legend() axes[2].grid(True, alpha0.3) plt.tight_layout() plt.savefig(microgrid_schedule.png, dpi150) plt.show() else: print(未找到最优解。请检查模型约束是否过紧如零碳约束导致不可行。)常见问题与排查模型不可行最常见原因是约束条件互相冲突。例如在零碳约束下不能买电如果某个时刻负荷远大于风光发电且储能电量不足功率平衡就无法满足。此时需要检查数据或者放宽零碳约束允许少量购电或引入碳交易。求解速度慢整数变量0-1变量会极大增加求解复杂度。如果调度周期T很长如96个15分钟间隔模型规模会变大。可以尝试a) 使用更高效的商业求解器Gurobi, CPLEXb) 松弛整数变量先求解连续问题看结果再引入整数约束c) 减少时间分辨率如从15分钟到1小时。结果不符合直觉比如储能一直在充放电“折腾”。可能是成本系数设置不合理或者缺少对储能循环寿命的惩罚。可以在目标函数中增加一项与充放电功率平方成正比的“磨损成本”来平滑储能动作。如何处理不确定性确定性模型假设预测是完美的。为了考虑风光和负荷的预测误差可以引入鲁棒优化。例如假设真实值在预测值±10%的区间内然后优化“最坏情况”下的成本。这会将问题转化为一个min-max问题求解更复杂但结果更稳健。另一种更实用的方法是滚动优化每执行完一个时段就用最新的实测数据更新预测重新求解后续时段实现闭环控制。5. 论文写作与模型拓展要点数模竞赛最终成果是论文。除了把模型和结果做出来如何清晰地呈现同样重要。对于C题医学图像论文摘要清晰说明问题、方法U-Net分割骨架化分析、关键指标分割Dice系数、狭窄识别准确率和结论。模型亮点强调数据预处理特别是CLAHE对医学图像的重要性说明损失函数选择DiceBCE如何解决类别不平衡讨论后处理中骨架化和直径测量的稳健性改进。结果展示务必提供可视化对比图包括原图、分割结果图、狭窄定位标注图。用表格列出在测试集上的量化评估指标。敏感性分析分析模型对关键超参数如学习率、损失函数权重的敏感性讨论狭窄判定阈值50%70%变化对分类结果的影响。对于D题微电网论文摘要阐明优化目标经济性、核心约束功率平衡、零碳、采用的方法混合整数线性规划和主要结论如零碳约束下成本增加XX%储能的最佳配置容量建议等。模型亮点详细阐述将物理问题转化为数学模型的每一步特别是整数变量充放电互斥的引入和零碳约束的数学表达。对比分析“有零碳约束”和“无零碳约束”两种场景下的调度策略和总成本突出零碳目标带来的挑战和储能价值的提升。结果分析除了展示调度计划图还要分析储能的“削峰填谷”作用计算风光渗透率、自平衡率等指标。用表格对比不同场景晴天/阴天、夏季/冬季下的优化结果。模型拓展讨论不确定性下的鲁棒优化模型或考虑设备投资成本储能容量、光伏装机的容量规划模型这可以作为“模型推广”部分体现思考的深度。给参赛者的最后建议数模竞赛是团队作战。C题和D题代表了两种典型风格。在组队时最好能涵盖编程能力强的队员实现算法、数学基础好的队员推导模型、写作表达能力强的队员撰写论文。拿到题目后花1-2小时进行头脑风暴和文献速查确定技术路线比盲目开始编程更重要。对于这类交叉学科题目快速学习新领域基础知识如冠脉造影临床知识、微电网基本架构的能力至关重要。最后永远记得用简洁的图表和严谨的数据来支撑你的论点这是打动评委的关键。 从医学影像分割到微电网优化：数模竞赛中的AI与运筹学实战
1. 项目概述当医学影像遇上能源优化看到“2026年河北省研究生数学建模 C/D 题”这个标题很多同学的第一反应可能是“又来画饼了”。确实预测两年后的赛题听起来有点玄乎但作为一名带过好几届数模竞赛、自己也从坑里爬出来的老手我得说这个标题背后透露出的命题趋势和备赛逻辑恰恰是现在准备竞赛最需要关注的核心。它把看似风马牛不相及的“冠脉造影狭窄识别”和“零碳园区微电网优化”放在一起绝不是随意拼凑而是精准地指向了当前交叉学科应用和复杂系统建模的两大前沿热点。简单来说这道“预测题”模拟了研究生数模竞赛的典型风格一道偏向数据科学与人工智能的C题一道偏向运筹优化与系统控制的D题。C题要求你从医学影像中自动识别血管狭窄这本质是一个计算机视觉中的图像分割与分类问题考验的是数据预处理、特征工程和模型构建能力。而D题则要求你对一个包含光伏、风电、储能和负荷的园区微电网进行优化调度目标是实现零碳运行与经济性的平衡这涉及到混合整数规划、动态优化甚至强化学习等运筹学方法。两者共同的核心都是将现实世界的复杂问题抽象为可计算的数学模型并给出可验证的解决方案。无论你是计算机、电气、自动化专业还是数学、统计背景这道题都能找到你的发力点。对于准备参赛的同学提前琢磨这类题目不是为了押题而是为了掌握一种“解题框架思维”——即面对一个陌生领域的具体问题如何快速进行问题分解、方法选型和模型构建。接下来我就结合自己多年踩坑和指导的经验对这两个题目进行深度拆解并提供一套可直接上手、能灵活调整的代码框架与建模思路。2. 核心思路拆解从问题到模型的映射逻辑面对一个综合性的数模题目最忌讳的就是一头扎进细节里开始编程。正确的打开方式是先进行顶层设计明确从“问题描述”到“最终论文”每一步要做什么用什么方法做以及为什么这么做。2.1 C题冠脉造影狭窄识别——一个标准的CV Pipeline医学影像分析题目在近年的国赛、美赛中屡见不鲜。冠脉造影狭窄识别的核心任务可以拆解为输入一张冠脉造影的X光影像输出图像中血管的轮廓以及狭窄部位的位置和程度评估。这本质上是一个“语义分割”“关键点/区域分类”的组合任务。为什么选择语义分割而不是单纯的目标检测因为血管是细长的管状结构狭窄表现为血管局部内径的变细。仅仅用检测框标出狭窄区域目标检测无法量化狭窄程度比如直径狭窄率。我们需要得到血管像素级的精确轮廓语义分割然后基于轮廓计算血管各处的直径从而找出直径小于正常值一定比例通常临床以50%或70%为界的狭窄段。整体建模Pipeline设计如下数据预处理与增强医学影像质量参差不齐存在对比度低、噪声大、造影剂不均匀等问题。预处理是关键第一步。血管分割模型采用编码器-解码器结构的神经网络如U-Net及其变体进行像素级分割将图像分为“血管”和“背景”两类。血管中心线提取与直径分析对分割出的二值血管掩膜图像采用骨架化算法提取单像素宽的中心线。沿着中心线法线方向测量血管的局部直径。狭窄识别与分类根据血管直径剖面识别局部最小值。计算狭窄率(1 - 狭窄处直径 / 参考正常直径) * 100%。根据狭窄率阈值如50%为轻度50%-70%为中度70%为重度对狭窄部位进行分类。结果可视化与评估将分割的血管轮廓、狭窄定位框和分类结果以热力图或标注框形式叠加在原图上并计算分割的Dice系数、狭窄识别的准确率/召回率等指标。注意临床实际中参考正常直径的选取有讲究。通常选取狭窄近端和远端看起来正常的血管段直径的平均值或采用同一血管的未狭窄段直径。在建模时如果数据未标注正常段可以假设血管直径变化平滑采用滑动窗口内的局部最大直径或拟合的基线直径作为参考。2.2 D题零碳园区微电网优化——一个典型的混合整数规划问题零碳园区微电网优化是一个多时间尺度、多设备耦合、多目标约束的优化问题。核心目标是在满足园区电负荷需求的前提下协调调度光伏发电、风力发电、蓄电池储能以及可能的外部电网购电假设外部电网为高碳电网购电会增加碳足迹使得一个调度周期如24小时内的总运行成本最低同时实现净零碳排放。为什么这是一个混合整数规划MIP问题因为模型中既包含连续变量如发电功率、储能充放电功率、荷电状态也包含整数变量或0-1变量。典型的整数变量场景包括储能的充放电状态同一时刻不能既充又放需要引入0-1变量表示模式、燃气轮机等可启停机组的启停状态、与电网的买卖电状态防止同时买卖等。这些“要么开要么关”的决策需要用整数变量来刻画。模型核心要素拆解决策变量各时刻光伏/风电出力通常视为不可控但可预测、储能充/放电功率、储能荷电状态SOC、从电网购电功率、向电网售电功率如果允许、可能的备用容量等。目标函数最小化总成本 购电成本 - 售电收入 设备运维成本与充放电量相关 储能折旧成本。为了实现“零碳”通常将碳排放量作为约束条件而非目标函数。约束条件功率平衡约束每一时刻发电总量 放电功率 购电功率 负荷 充电功率 售电功率。设备运行约束光伏/风电出力上下限储能SOC上下限、充放电功率上下限、充放电效率、初始/终止SOC约束充放电互斥约束需0-1变量。电网交互约束购电/售电功率上下限可能存在的功率因数约束。零碳约束调度周期内总碳排放量 ≤ 0。碳排放来源主要是外购电力假设其有碳强度系数如kgCO₂/kWh。园区内光伏、风电、储能视为零碳。因此约束可转化为外购电量 * 碳强度系数 ≤ 0这实质上强制要求外购电量为零除非题目允许碳交易即购买绿证或碳汇来抵消。这是一个非常强的约束会极大改变优化策略。不确定性处理光伏、风电和负荷预测都存在误差。高级的模型会考虑鲁棒优化或随机规划引入不确定性集合或场景使调度方案对预测误差不敏感。对于初次接触可以先做确定性优化。3. C题代码框架与实操要点这里提供一个基于Python和PyTorch的简化版代码框架重点展示核心流程和易错点。3.1 数据准备与预处理模块医学影像数据通常以DICOM格式存储但竞赛为了方便很可能提供的是处理后的JPEG或PNG图像甚至可能是公开数据集如部分标注的冠脉造影数据集。import cv2 import numpy as np import torch from torch.utils.data import Dataset, DataLoader import albumentations as A from albumentations.pytorch import ToTensorV2 class CoronaryAngioDataset(Dataset): def __init__(self, image_paths, mask_pathsNone, is_trainTrue, img_size512): self.image_paths image_paths self.mask_paths mask_paths self.is_train is_train self.img_size img_size # 训练和验证/测试采用不同的数据增强策略 if self.is_train: self.transform A.Compose([ A.Resize(heightimg_size, widthimg_size), A.RandomRotate90(p0.5), A.HorizontalFlip(p0.5), A.VerticalFlip(p0.5), # 医学影像增强对比度受限自适应直方图均衡化(CLAHE)非常有效 A.CLAHE(clip_limit2.0, tile_grid_size(8,8), p0.8), A.RandomBrightnessContrast(brightness_limit0.1, contrast_limit0.1, p0.5), A.GaussNoise(var_limit(10.0, 50.0), p0.3), A.Normalize(mean[0.485], std[0.229]), # 灰度图单通道 ToTensorV2(), ]) else: self.transform A.Compose([ A.Resize(heightimg_size, widthimg_size), A.CLAHE(clip_limit2.0, tile_grid_size(8,8), always_applyTrue), # 测试时也做CLAHE A.Normalize(mean[0.485], std[0.229]), ToTensorV2(), ]) def __getitem__(self, idx): image cv2.imread(self.image_paths[idx], cv2.IMREAD_GRAYSCALE) if image is None: raise FileNotFoundError(f无法读取图像: {self.image_paths[idx]}) # 确保图像是二维的 if len(image.shape) 3: image cv2.cvtColor(image, cv2.COLOR_BGR2GRAY) if self.mask_paths is not None: mask cv2.imread(self.mask_paths[idx], cv2.IMREAD_GRAYSCALE) mask (mask 127).astype(np.uint8) # 二值化 augmented self.transform(imageimage, maskmask) image_tensor augmented[image] mask_tensor augmented[mask].float().unsqueeze(0) # 增加通道维 return image_tensor, mask_tensor else: augmented self.transform(imageimage) return augmented[image], self.image_paths[idx] # 测试时返回路径用于保存结果 def __len__(self): return len(self.image_paths)实操心得CLAHE是关键对于对比度低的X光影像CLAHE比简单的直方图均衡化效果更好能增强血管与背景的对比且不易过度增强噪声。谨慎使用几何变换旋转、翻转对于血管结构是有效的因为血管方向不固定。但像弹性形变这类强变换要小心可能扭曲血管的拓扑结构如分叉点。二值化阈值标注的掩膜可能是0-255的灰度图需要统一二值化。127是个常用阈值但最好检查一下标注数据的实际范围。3.2 模型构建U-Net及其改进U-Net是医学分割的经典网络。这里实现一个基础版并讨论几个改进方向。import torch.nn as nn import torch.nn.functional as F class DoubleConv(nn.Module): (卷积 BN ReLU) * 2 def __init__(self, in_channels, out_channels): super().__init__() self.double_conv nn.Sequential( nn.Conv2d(in_channels, out_channels, kernel_size3, padding1), nn.BatchNorm2d(out_channels), nn.ReLU(inplaceTrue), nn.Conv2d(out_channels, out_channels, kernel_size3, padding1), nn.BatchNorm2d(out_channels), nn.ReLU(inplaceTrue) ) def forward(self, x): return self.double_conv(x) class UNet(nn.Module): def __init__(self, n_channels1, n_classes1): super(UNet, self).__init__() self.n_channels n_channels self.n_classes n_classes self.inc DoubleConv(n_channels, 64) self.down1 nn.Sequential(nn.MaxPool2d(2), DoubleConv(64, 128)) self.down2 nn.Sequential(nn.MaxPool2d(2), DoubleConv(128, 256)) self.down3 nn.Sequential(nn.MaxPool2d(2), DoubleConv(256, 512)) self.down4 nn.Sequential(nn.MaxPool2d(2), DoubleConv(512, 1024)) self.up1 nn.ConvTranspose2d(1024, 512, kernel_size2, stride2) self.conv1 DoubleConv(1024, 512) # 1024 512(skip) 512(up) self.up2 nn.ConvTranspose2d(512, 256, kernel_size2, stride2) self.conv2 DoubleConv(512, 256) self.up3 nn.ConvTranspose2d(256, 128, kernel_size2, stride2) self.conv3 DoubleConv(256, 128) self.up4 nn.ConvTranspose2d(128, 64, kernel_size2, stride2) self.conv4 DoubleConv(128, 64) self.outc nn.Conv2d(64, n_classes, kernel_size1) def forward(self, x): x1 self.inc(x) x2 self.down1(x1) x3 self.down2(x2) x4 self.down3(x3) x5 self.down4(x4) x self.up1(x5) # 拼接跳跃连接注意空间尺寸要对齐 x torch.cat([x, x4], dim1) x self.conv1(x) x self.up2(x) x torch.cat([x, x3], dim1) x self.conv2(x) x self.up3(x) x torch.cat([x, x2], dim1) x self.conv3(x) x self.up4(x) x torch.cat([x, x1], dim1) x self.conv4(x) logits self.outc(x) return logits # 改进方向讨论 # 1. 深度监督在中间层也添加辅助输出帮助梯度回传缓解梯度消失。 # 2. 注意力机制在跳跃连接处加入注意力门Attention Gate让网络更关注血管区域抑制背景噪声。 # 3. 使用预训练编码器将U-Net的编码器替换为ResNet、EfficientNet等在大规模图像数据上预训练的网络可以提升特征提取能力尤其在小数据集上效果显著。 # 4. 损失函数二值分割常用Dice Loss BCE Loss的组合。Dice Loss直接优化分割区域的重叠度对类别不平衡血管像素少问题更鲁棒。3.3 后处理与狭窄分析模块这是将分割结果转化为最终答案的关键也是容易出bug的地方。import cv2 import numpy as np from skimage import morphology, measure def analyze_stenosis_from_mask(binary_mask): 从二值血管掩膜分析狭窄情况。 参数: binary_mask: numpy数组二值图像1为血管0为背景。 返回: results: 列表每个元素是一个字典描述一个狭窄区域。 # 1. 细化骨架化获取中心线 # 使用中轴骨架化保留拓扑结构 skeleton morphology.skeletonize(binary_mask 0) # 2. 获取中心线的坐标点 skeleton_points np.column_stack(np.where(skeleton 0)) # 3. 构建中心线的图结构用于追踪血管分支 # 这里简化处理假设血管是单连通的无分叉或处理主干。 # 对于复杂分叉需要先对骨架进行端点、交叉点分析分割成多个段。 from scipy.spatial import KDTree tree KDTree(skeleton_points) # 4. 沿着中心线测量局部直径 diameters [] for i, point in enumerate(skeleton_points): # 对于中心线上每个点计算到最近背景点的距离*2 作为该点处血管直径的近似 # 更准确的做法在该点法线方向需估算骨架方向上向两边搜索背景点。 # 这里使用距离变换作为简化方法。 pass # 具体实现需要结合距离变换 # 5. 平滑直径曲线并寻找局部极小值狭窄点 from scipy.signal import savgol_filter, find_peaks smoothed_diameters savgol_filter(diameters, window_length11, polyorder3) # 寻找局部最小值狭窄注意find_peaks找的是极大值对负值取反找极小值 minima_indices, _ find_peaks(-smoothed_diameters, prominence0.5) # prominence需要根据数据调整 # 6. 计算狭窄率并分类 results [] for idx in minima_indices: stenosis_diameter smoothed_diameters[idx] # 参考直径取狭窄点前后一段正常血管直径的平均值例如前后各10个点排除其他狭窄点 window_size 10 start max(0, idx - window_size) end min(len(smoothed_diameters), idx window_size 1) # 简单起见取窗口内最大值作为参考直径假设血管直径变化平缓 reference_diameter np.max(smoothed_diameters[start:end]) stenosis_ratio (1 - stenosis_diameter / reference_diameter) * 100 severity 轻度 if stenosis_ratio 70: severity 重度 elif stenosis_ratio 50: severity 中度 # 获取狭窄点在原图中的坐标近似 stenosis_point skeleton_points[idx] results.append({ location_xy: (int(stenosis_point[1]), int(stenosis_point[0])), # (x, y) stenosis_ratio_percent: round(stenosis_ratio, 2), severity: severity, reference_diameter: round(reference_diameter, 2), stenosis_diameter: round(stenosis_diameter, 2) }) return results # 可视化函数 def visualize_results(original_image, binary_mask, stenosis_results, save_pathNone): vis_img cv2.cvtColor(original_image, cv2.COLOR_GRAY2BGR) # 用半透明红色覆盖血管区域 mask_color np.zeros_like(vis_img) mask_color[binary_mask 0] [0, 0, 255] # BGR中的红色 cv2.addWeighted(mask_color, 0.3, vis_img, 0.7, 0, vis_img) # 标记狭窄点 for res in stenosis_results: x, y res[location_xy] cv2.circle(vis_img, (x, y), radius5, color(0, 255, 0), thickness-1) # 绿色圆点 label f{res[severity]}({res[stenosis_ratio_percent]}%) cv2.putText(vis_img, label, (x10, y-10), cv2.FONT_HERSHEY_SIMPLEX, 0.5, (0, 255, 0), 1) if save_path: cv2.imwrite(save_path, vis_img) return vis_img踩坑实录骨架化噪声二值掩膜边缘如果有毛刺骨架化会产生大量短小分支干扰中心线提取。务必先对掩膜进行morphology.opening开运算去除小噪声。血管分叉处理真实的冠脉血管是有分叉的。简单的骨架化会在分叉点产生一个“星形”结构导致中心线追踪混乱。需要在骨架化后进行端点、交叉点检测将血管分解成多个独立的“段”再分别分析。可以使用skimage.morphology中的remove_small_objects和skan库Skeleton Analysis进行更专业的骨架分析。直径测量不准在弯曲的血管处法线方向估算不准会导致直径测量偏差。一个更稳健的方法是使用距离变换dist_transform cv2.distanceTransform(binary_mask.astype(np.uint8), cv2.DIST_L2, 3)。对于中心线上的点其距离变换值的两倍可以较好地近似该点处的血管半径。但要注意在血管边缘不平滑时仍有误差。狭窄判定阈值临床上的50%、70%阈值是金标准但在模型评估时如果数据标注的“狭窄”是基于视觉评估的可能与计算的直径狭窄率有系统偏差。建议在验证集上做一个简单的校准看看模型计算的狭窄率与标注的狭窄等级如果有的话之间的对应关系微调阈值。4. D题建模与求解框架微电网优化建模通常使用线性规划LP或混合整数线性规划MILP工具。Python中PuLP和CVXPY对于线性问题是易上手的建模库而Gurobi、CPLEX是强大的商业求解器学术可免费申请license。这里以PuLP为例构建一个确定性优化模型。4.1 模型参数与变量定义假设调度周期T24小时时间间隔Δt1小时。import pulp import numpy as np # 参数定义 (示例值需根据题目数据修改) T 24 # 时间步数 # 预测数据 (kW) P_pv_pred np.array([...]) # 光伏预测出力长度T P_wind_pred np.array([...]) # 风电预测出力 P_load_pred np.array([...]) # 电负荷预测 # 电价 (元/kWh) price_buy np.array([...]) # 购电价分时电价 price_sell np.array([...]) # 售电价通常低于购电价 carbon_intensity_grid 0.8 # 电网碳排放强度 kgCO2/kWh # 储能参数 E_batt_max 500 # 储能最大容量 kWh E_batt_min 50 # 储能最小容量防止过放设为最大容量的10% P_charge_max 100 # 最大充电功率 kW P_discharge_max 100 # 最大放电功率 kW eta_charge 0.95 # 充电效率 eta_discharge 0.95 # 放电效率 E_batt_initial 250 # 初始电量 kWh E_batt_end 250 # 调度周期末期望电量 kWh通常等于初始值 # 创建问题 prob pulp.LpProblem(ZeroCarbon_Microgrid_Optimization, pulp.LpMinimize) # 决策变量 # 连续变量 P_grid_buy pulp.LpVariable.dicts(P_grid_buy, range(T), lowBound0) # 购电功率 P_grid_sell pulp.LpVariable.dicts(P_grid_sell, range(T), lowBound0) # 售电功率 P_batt_charge pulp.LpVariable.dicts(P_batt_charge, range(T), lowBound0, upBoundP_charge_max) P_batt_discharge pulp.LpVariable.dicts(P_batt_discharge, range(T), lowBound0, upBoundP_discharge_max) E_batt pulp.LpVariable.dicts(E_batt, range(T1), lowBoundE_batt_min, upBoundE_batt_max) # 从t0到tT # 0-1整数变量 (确保不同时充放电) u_charge pulp.LpVariable.dicts(u_charge, range(T), catBinary) # 1表示充电 u_discharge pulp.LpVariable.dicts(u_discharge, range(T), catBinary) # 1表示放电 # 目标函数最小化总运行成本 # 成本 购电成本 - 售电收入 储能循环损耗成本简化为与充放电量成正比 cost_buy pulp.lpSum([P_grid_buy[t] * price_buy[t] for t in range(T)]) cost_sell pulp.lpSum([P_grid_sell[t] * price_sell[t] for t in range(T)]) # 假设储能充放电损耗成本系数为0.05元/kWh cost_batt_cycle 0.05 * pulp.lpSum([P_batt_charge[t] P_batt_discharge[t] for t in range(T)]) prob cost_buy - cost_sell cost_batt_cycle4.2 约束条件构建这是模型的核心需要严谨地表达所有物理和运行限制。# 约束条件 for t in range(T): # 1. 功率平衡约束 (核心) prob (P_pv_pred[t] P_wind_pred[t] P_batt_discharge[t] P_grid_buy[t] P_load_pred[t] P_batt_charge[t] P_grid_sell[t]), fPower_Balance_t{t} # 2. 储能电量动态方程 if t 0: prob (E_batt[t] E_batt_initial), fBatt_Init prob (E_batt[t1] E_batt[t] P_batt_charge[t] * eta_charge * 1 # Δt1小时能量功率*时间 - P_batt_discharge[t] / eta_discharge * 1), fBatt_Dynamics_t{t} # 3. 充放电互斥约束 (Big-M法) M max(P_charge_max, P_discharge_max) * 1.5 # 一个足够大的数 prob (P_batt_charge[t] u_charge[t] * M), fCharge_Logic_1_t{t} prob (P_batt_discharge[t] u_discharge[t] * M), fDischarge_Logic_1_t{t} prob (u_charge[t] u_discharge[t] 1), fCharge_Discharge_Mutex_t{t} # 4. 电网交互约束 (假设最大交互功率为200kW) P_grid_max 200 prob (P_grid_buy[t] P_grid_max), fGrid_Buy_Max_t{t} prob (P_grid_sell[t] P_grid_max), fGrid_Sell_Max_t{t} # 通常不允许同时购售电但上述互斥约束已通过功率平衡和成本函数自然实现因为售电价格低同时购售电不经济。如需严格约束可再加0-1变量。 # 储能周期首尾电量约束保证可持续运行 prob (E_batt[T] E_batt_end), Batt_End_Min # 也可设为等于强制要求周期末回到初始状态 # prob (E_batt[T] E_batt_end), Batt_End_Equal # 零碳约束 # 总碳排放 总购电量 * 碳强度 total_carbon_emission pulp.lpSum([P_grid_buy[t] for t in range(T)]) * carbon_intensity_grid # 强制零碳总碳排放 0 prob (total_carbon_emission 0), Zero_Carbon_Constraint # 注意这个约束非常强意味着 P_grid_buy[t] 对所有t都必须为0。 # 如果题目允许少量碳排放或碳交易约束应改为 Carbon_Budget (碳排放预算)。4.3 模型求解与结果分析# 求解 # 指定求解器例如CBC开源或Gurobi需安装且配置license solver pulp.PULP_CBC_CMD(msgFalse) # 使用CBC不输出求解日志 # solver pulp.GUROBI_CMD(msgTrue) # 如果安装了Gurobi prob.solve(solver) # 输出结果 print(f求解状态: {pulp.LpStatus[prob.status]}) print(f最优总成本: {pulp.value(prob.objective):.2f} 元) if prob.status pulp.LpOptimal: # 提取变量值 P_buy_opt [pulp.value(P_grid_buy[t]) for t in range(T)] P_sell_opt [pulp.value(P_grid_sell[t]) for t in range(T)] P_charge_opt [pulp.value(P_batt_charge[t]) for t in range(T)] P_discharge_opt [pulp.value(P_batt_discharge[t]) for t in range(T)] E_batt_opt [pulp.value(E_batt[t]) for t in range(T1)] # 计算实际碳排放 total_buy sum(P_buy_opt) actual_carbon total_buy * carbon_intensity_grid print(f总购电量: {total_buy:.2f} kWh, 理论碳排放: {actual_carbon:.2f} kgCO2) # 在零碳约束下total_buy应为0 # 可视化调度计划 import matplotlib.pyplot as plt fig, axes plt.subplots(3, 1, figsize(12, 10)) t_range list(range(T)) axes[0].plot(t_range, P_load_pred, label负荷, colorblack, linewidth2) axes[0].plot(t_range, P_pv_pred, label光伏, colorgold) axes[0].plot(t_range, P_wind_pred, label风电, colorcyan) axes[0].set_ylabel(功率 (kW)) axes[0].set_title(预测的发电与负荷) axes[0].legend() axes[0].grid(True, alpha0.3) axes[1].bar(t_range, P_charge_opt, label充电, colorgreen, alpha0.6) axes[1].bar(t_range, [-x for x in P_discharge_opt], label放电, colorred, alpha0.6) # 放电为负 axes[1].plot(t_range, E_batt_opt[:-1], label储能电量, colorblue, markero) axes[1].set_ylabel(功率/电量 (kW/kWh)) axes[1].set_title(储能调度) axes[1].legend() axes[1].grid(True, alpha0.3) axes[2].bar(t_range, P_buy_opt, label购电, colororange, alpha0.6) axes[2].bar(t_range, [-x for x in P_sell_opt], label售电, colorpurple, alpha0.6) # 售电为负 axes[2].set_xlabel(时间 (h)) axes[2].set_ylabel(功率 (kW)) axes[2].set_title(与电网交互) axes[2].legend() axes[2].grid(True, alpha0.3) plt.tight_layout() plt.savefig(microgrid_schedule.png, dpi150) plt.show() else: print(未找到最优解。请检查模型约束是否过紧如零碳约束导致不可行。)常见问题与排查模型不可行最常见原因是约束条件互相冲突。例如在零碳约束下不能买电如果某个时刻负荷远大于风光发电且储能电量不足功率平衡就无法满足。此时需要检查数据或者放宽零碳约束允许少量购电或引入碳交易。求解速度慢整数变量0-1变量会极大增加求解复杂度。如果调度周期T很长如96个15分钟间隔模型规模会变大。可以尝试a) 使用更高效的商业求解器Gurobi, CPLEXb) 松弛整数变量先求解连续问题看结果再引入整数约束c) 减少时间分辨率如从15分钟到1小时。结果不符合直觉比如储能一直在充放电“折腾”。可能是成本系数设置不合理或者缺少对储能循环寿命的惩罚。可以在目标函数中增加一项与充放电功率平方成正比的“磨损成本”来平滑储能动作。如何处理不确定性确定性模型假设预测是完美的。为了考虑风光和负荷的预测误差可以引入鲁棒优化。例如假设真实值在预测值±10%的区间内然后优化“最坏情况”下的成本。这会将问题转化为一个min-max问题求解更复杂但结果更稳健。另一种更实用的方法是滚动优化每执行完一个时段就用最新的实测数据更新预测重新求解后续时段实现闭环控制。5. 论文写作与模型拓展要点数模竞赛最终成果是论文。除了把模型和结果做出来如何清晰地呈现同样重要。对于C题医学图像论文摘要清晰说明问题、方法U-Net分割骨架化分析、关键指标分割Dice系数、狭窄识别准确率和结论。模型亮点强调数据预处理特别是CLAHE对医学图像的重要性说明损失函数选择DiceBCE如何解决类别不平衡讨论后处理中骨架化和直径测量的稳健性改进。结果展示务必提供可视化对比图包括原图、分割结果图、狭窄定位标注图。用表格列出在测试集上的量化评估指标。敏感性分析分析模型对关键超参数如学习率、损失函数权重的敏感性讨论狭窄判定阈值50%70%变化对分类结果的影响。对于D题微电网论文摘要阐明优化目标经济性、核心约束功率平衡、零碳、采用的方法混合整数线性规划和主要结论如零碳约束下成本增加XX%储能的最佳配置容量建议等。模型亮点详细阐述将物理问题转化为数学模型的每一步特别是整数变量充放电互斥的引入和零碳约束的数学表达。对比分析“有零碳约束”和“无零碳约束”两种场景下的调度策略和总成本突出零碳目标带来的挑战和储能价值的提升。结果分析除了展示调度计划图还要分析储能的“削峰填谷”作用计算风光渗透率、自平衡率等指标。用表格对比不同场景晴天/阴天、夏季/冬季下的优化结果。模型拓展讨论不确定性下的鲁棒优化模型或考虑设备投资成本储能容量、光伏装机的容量规划模型这可以作为“模型推广”部分体现思考的深度。给参赛者的最后建议数模竞赛是团队作战。C题和D题代表了两种典型风格。在组队时最好能涵盖编程能力强的队员实现算法、数学基础好的队员推导模型、写作表达能力强的队员撰写论文。拿到题目后花1-2小时进行头脑风暴和文献速查确定技术路线比盲目开始编程更重要。对于这类交叉学科题目快速学习新领域基础知识如冠脉造影临床知识、微电网基本架构的能力至关重要。最后永远记得用简洁的图表和严谨的数据来支撑你的论点这是打动评委的关键。