简介面向需要求解多目标优化问题的 MATLAB 用户这份资源将遗传算法实现、标准测试函数与工具文档整合在一起帮助读者理解帕累托最优解集的核心概念并快速搭建自己的优化流程。资源包共 64 个文件其中 60 个为脚本或函数文件覆盖种群初始化、适应度排序、选择、交叉与变异等关键算子并附有主程序、工具箱说明文档、文本说明及一个 PostScript 格式图片整体压缩后仅 2.74 MB。借助内置的多个测试函数和可视化工具配合遗传算法工具箱里丰富的算子策略可以先运行主程序熟悉完整求解流程再替换成自己的目标函数并设置约束观察帕累托前沿的收敛过程非常适用于数学建模竞赛、课程设计和工程多方案比选。整套代码结构清晰函数拆分合理便于按模块理解与二次开发。目前已有 2540 人学习浏览值得中高级 MATLAB 用户仔细研读算法细节、调整参数从而获得更优的妥协解。1. 多目标优化落地用 MATLAB 跑出 Pareto 前沿的那点事帮某导师做设备参数寻优时遇到过一件反直觉的事把所有目标加权成一个数在单目标里找最优跑了一周也没给出一组让所有人满意的参数改用多目标优化之后一上午拿回一组 Pareto 解各目标互有取舍剩下的事交给了拍板的人。这就是这套多目标函数优化 MATLAB 工具的价值——它不追求“唯一最优解”而是用非支配排序和拥挤度机制把一组互不占优的候选解摊开在一个平面上让你用眼睛做决策。适合拿来做算法对比、课题起步或者给真实工程问题换一种寻优姿势。读完这篇你至少能做到两件事把源码跑起来画出一张 Pareto 前沿图以及把测试函数换成你自己的目标函数。2. 架构先行从种群初始化到 Pareto 输出的模块划分2.1 为什么选 MATLAB 而不是 Python做多目标优化Python 生态里不是没有 pymoo、DEAP 这类库但 MATLAB 在这个场景里有三个不可替代的体验第一画图即所得plot一敲 Pareto 前沿就在屏幕上交互式旋转三维散点图也顺手对“先看图、再调参、再判断”的优化流程极其友好第二矩阵运算是语言内建能力种群在 MATLAB 里天然是矩阵向量化写法和数值计算的老底子都在第三学术语境里 MATLAB 出图方便论文图例、坐标轴标注、字体控制成熟做完实验直接导出成插图。代价是速度。for 循环一多MATLAB 立刻变“慢缸”。所以这套代码里凡是能向量化的地方非支配判定的批量比较、拥挤度排序都不该写三层嵌套。资源里如果把纯循环版本和向量化版本都放进来优先看后者——性能差距在 500 个体、500 代时就非常可观前者可能跑到怀疑人生。2.2 目录结构与文件职责这类资源拿到手我一般先不看 README直接tree扫一遍结构。标准的多目标优化 MATLAB 项目通常按“算法核心 / 测试问题 / 实验脚本”三层组织和你看到的源码包基本一致multi-objective-optimization-matlab/ ├── src/ # 算法核心 │ ├── nsga2.m # 主循环遗传迭代与精英保留 │ ├── fastNonDominatedSort.m # 非支配排序输出分层 rank │ ├── crowdingDistance.m # 拥挤度计算维持多样性 │ ├── tournamentSelect.m # 二元锦标赛选择 │ ├── sbxCrossover.m # 模拟二进制交叉 SBX │ ├── polyMutation.m # 多项式变异 │ └── environmentalSelection.m # 环境选择合并后择优 ├── problems/ # 目标函数 │ ├── zdt1.m # ZDT1凸 Pareto 前沿 │ ├── zdt2.m # ZDT2非凸前沿 │ └── zdt3.m # ZDT3不连续前沿 ├── scripts/ │ ├── run_zdt1.m # 一键跑通实验 │ └── compare_zdt.m # 不同代数结果对比 └── README.md逐文件说nsga2.m是整个包的骨架它负责把初始化、评估、排序、选择、交叉、变异、环境选择串起来大部分时候你不需要改它只需要知道它调用了哪些函数。fastNonDominatedSort.m和crowdingDistance.m是算法的心脏前者管“谁更优”后者管“别挤在一起”。problems/目录是你二次开发时主要动手的地方——把 ZDT1 换成你自己的目标函数项目就能从教学示例变成工程工具。scripts/里的脚本是入口跑实验从这些文件进。2.3 参数约定与默认值多目标优化最怕“一套参数走天下”。NSGA-II 的核心参数集中在nsga2.m顶部的结构体里常见的默认配置如下参数符号默认值影响种群规模popSize100太小时 Pareto 前沿稀疏太大时单代耗时上升迭代代数gen200控制搜索深度判断是否收敛要看前沿变化交叉概率pc0.9产生子代的主要手段变异概率pm1 / dim每个决策变量的变异概率通常取维度倒数交叉分布指数etaC20越大子代越接近父代搜索越“保守”变异分布指数etaM20控制变异扰动的步长范围决策变量边界lb / ubzeros / ones必须与测试问题的定义域一致常被忽略我一般会把lb、ub从问题函数里提出来而不是写死在主循环里。这样换测试问题或接真实工程变量时只需改一处。记住一个原则改任何参数前先确认当前问题是什么形状的前沿再决定要不要动种群规模——这个后面第 5 章会展开说。3. 核心算法拆解非支配排序与拥挤度如何撑起 NSGA-II3.1 种群初始化与目标评估NSGA-II 从随机种群起步初始化时保证每个决策变量均匀落在定义域内function pop initializePopulation(popSize, dim, lb, ub) % 生成 popSize x dim 的种群每行是一个个体 pop rand(popSize, dim) .* (ub - lb) lb; endrand生成 (0,1) 均匀分布随机数乘上(ub - lb)再偏移lb每个变量都在边界内。这一步没有什么玄学关键是维度dim必须和测试函数的输入长度一致。ZDT1 常用 30 维决策变量如果你把dim写成 10前沿形状会不一样但也能跑容易误判算法效果。目标评估在项目里通常写成独立函数输入一个个体行向量输出一行目标值function f zdt1(x) % ZDT1 测试函数双目标最小化 n length(x); g 1 9 * sum(x(2:n)) / (n - 1); % 收敛性函数 f1 x(1); % 第一个目标只依赖第一个变量 f2 g * (1 - sqrt(f1 / g)); % 第二个目标与所有变量相关 f [f1, f2]; % 返回行向量双目标 end这里有一个接口约定值得注意目标函数返回的是行向量不是列向量也不是矩阵。后面fastNonDominatedSort会把整个种群的目标值堆成 N×M 矩阵每一行是一个个体每一列是一个目标。如果目标函数返回了列向量排序结果会完全错乱这是新手最容易踩的坑之一。3.2 非支配排序把种群分层非支配排序是 NSGA-II 的“裁判”。对两个个体 A、B若 A 的所有目标都不差于 B且至少一个目标严格优于 B则称 A 支配 B。排序的目标是把种群分成若干层第一层是没被任何人支配的个体第二层是除去第一层后没被支配的以此类推。经典实现如下function [rank, F] fastNonDominatedSort(pop) % pop: N x M 目标矩阵M 2 N size(pop, 1); n zeros(N, 1); % n(i) 表示支配个体 i 的解的数量 S cell(N, 1); % S{i} 存储被个体 i 支配的解的下标 for i 1:N for j 1:N if i j, continue; end % 判断 i 是否支配 j if all(pop(i,:) pop(j,:)) any(pop(i,:) pop(j,:)) S{i} [S{i}, j]; % 判断 j 是否支配 i elseif all(pop(j,:) pop(i,:)) any(pop(j,:) pop(i,:)) n(i) n(i) 1; end end end F{1} find(n 0); % 第一层不被任何个体支配 k 1; while ~isempty(F{k}) Q []; % 暂存下一层 for i F{k} for j S{i} n(j) n(j) - 1; if n(j) 0, Q [Q, j]; end end end k k 1; if isempty(Q), break; end F{k} Q; end rank zeros(N, 1); for k 1:length(F) rank(F{k}) k; end end三个关键细节all(pop(i,:) pop(j,:))配合any(pop(i,:) pop(j,:))排除了“完全相同”的个体互相支配的情况——两个目标值完全一样的个体算互相不支配落在同一层这是合理行为。S用 cell 数组存下被支配集合后续减计数时直接遍历避免重复比较。最外层的双重循环复杂度是 O(MN²)种群 500 以内还能忍超过 1000 就明显卡顿优化办法是先用按某个目标排序的剪枝策略替代双层循环后面避坑章再细说。3.3 拥挤度计算让解铺开而不是挤成一团非支配排序保证“收敛”但光收敛不够——如果所有解都堆在 Pareto 前沿的一个角落这个结果对决策者没有意义。拥挤度就是惩罚“扎堆”的指标每个个体的拥挤度等于它周围解的密度越大代表越“孤独”在环境选择中越有优势function dist crowdingDistance(pop) % pop: N x M 目标矩阵 N size(pop, 1); M size(pop, 2); dist zeros(N, 1); for m 1:M [~, idx] sort(pop(:, m)); % 按第 m 个目标排序 dist(idx(1)) Inf; % 边界点永远保留 dist(idx(end)) Inf; fmin pop(idx(1), m); fmax pop(idx(end), m); if fmax fmin, continue; end % 该目标所有值相同跳过 for i 2:N-1 dist(idx(i)) dist(idx(i)) ... (pop(idx(i1), m) - pop(idx(i-1), m)) / (fmax - fmin); end end end逐行解释先按某个目标升序排列两端的边界解拥挤度直接设为Inf保证他们永远优先入选这保持了 Pareto 前沿的“跨度”。内部解的拥挤度用“前后邻居在该目标上的差除以该目标的总范围”来归一化最后把每个目标上的贡献累加。注意fmax fmin的短路判断如果某个目标所有个体取值相同比如收敛到一点直接跳过否则除零会得到 NaN污染整个排序结果。拥挤度的物理含义可以这么理解它度量以该个体为中心的长方形周长。前沿两端是“无限宽”中间解越稀疏周长越大淘汰时先淘汰拥挤度小的——也就是被包围得最严实的那些解。3.4 锦标赛选择与遗传算子选择、交叉、变异是 NSGA-II 的“变异引擎”。选择采用二元锦标赛随机抽两个个体优先选 rank 小的rank 相同时选拥挤度大的function idx tournamentSelect(rank, dist, k) % 从种群中随机抽 k 个候选返回最优者下标 candidates randperm(length(rank), k); bestPos 1; for i 2:k cr rank(candidates(i)); dr dist(candidates(i)); cb rank(candidates(bestPos)); db dist(candidates(bestPos)); if cr cb || (cr cb dr db) bestPos i; end end idx candidates(bestPos); endk默认取 2称为二元锦标赛。两个候选不重复randperm保证这一点。这个函数只返回一个父本下标主循环里调用两次得到两个不同的父本做交叉。交叉用模拟二进制交叉 SBX它模拟二进制编码的交叉效果但直接作用在实数上function [c1, c2] sbxCrossover(p1, p2, pc, etaC, lb, ub) c1 p1; c2 p2; if rand pc, return; end % 未命中交叉概率子代等于父代 for i 1:length(p1) if rand 0.5 % 每个变量按概率决定是否交叉 if abs(p1(i) - p2(i)) 1e-10 continue; % 两父本该变量几乎相同跳过 end u rand; if u 0.5 beta (2*u)^(1/(etaC1)); else beta (1/(2*(1-u)))^(1/(etaC1)); end c1(i) 0.5*((1beta)*p1(i) (1-beta)*p2(i)); c2(i) 0.5*((1-beta)*p1(i) (1beta)*p2(i)); % 越界裁剪确保子代仍在定义域内 c1(i) min(max(c1(i), lb(i)), ub(i)); c2(i) min(max(c2(i), lb(i)), ub(i)); end end endetaC是交叉分布指数值越大beta越接近 1子代越像父本搜索越局部值越小子代越可能远离父本探索性越强。工程上取 20 是原论文推荐值改了它等于改了算法的探索-开发平衡。多项式变异同理etaM控制扰动幅度。这两个算子本身不复杂复杂的是它们的参数要和种群规模、代数配合起来调单拎出来调其中一个没有意义。3.5 环境选择与主循环NSGA-II 的精髓在“精英保留”每一代把父代和子代合并成 2N 规模的大池子按 (rank, 拥挤度) 排序选前 N 个进入下一代。这样做的好处是优秀个体永远不会丢单调不退化function nextPop environmentalSelection(pop, rank, dist, popSize) nextPop []; count 0; for r 1:max(rank) idx find(rank r); if count length(idx) popSize % 整层都能装下直接全部保留 nextPop [nextPop; pop(idx, :)]; count count length(idx); else % 这一层只能装一部分按拥挤度降序取头部 [~, order] sort(dist(idx), descend); need popSize - count; nextPop [nextPop; pop(idx(order(1:need)), :)]; break; end end end主循环把这些零件串起来function [pop, obj, rank] nsga2(problem, params) pop initializePopulation(params.popSize, params.dim, params.lb, params.ub); for g 1:params.gen obj evaluateAll(problem, pop); % 批量计算目标值 [rank, ~] fastNonDominatedSort(obj); dist crowdingDistance(obj); % 生成子代 offspring []; while size(offspring, 1) params.popSize p1 tournamentSelect(rank, dist, 2); p2 tournamentSelect(rank, dist, 2); [c1, c2] sbxCrossover(pop(p1,:), pop(p2,:), params.pc, params.etaC, params.lb, params.ub); c1 polyMutation(c1, params.pm, params.etaM, params.lb, params.ub); c2 polyMutation(c2, params.pm, params.etaM, params.lb, params.ub); offspring [offspring; c1; c2]; end offspring offspring(1:params.popSize, :); % 去掉多出来的一个若为奇数 % 合并父代与子代环境选择 allPop [pop; offspring]; allObj evaluateAll(problem, allPop); [rank, ~] fastNonDominatedSort(allObj); dist crowdingDistance(allObj); pop environmentalSelection(allPop, rank, dist, params.popSize); end obj evaluateAll(problem, pop); end循环结构不复杂但有三处容易写错子代 while 循环每次生成 2 个个体popSize为奇数时会多出一个所以有offspring(1:popSize,:)这行兜底环境选择输入的是合并后的大种群rank 和拥挤度必须重新算不能复用子代生成前的旧值每一代结束更新pop后下一轮迭代开头要先重新算目标值因为pop的内容已经变了。4. 跑通实验用 ZDT1 验证实现再和理论前沿对齐4.1 运行脚本20 行代码获得第一张 Pareto 图资源里的scripts/run_zdt1.m就是入口核心逻辑如下% run_zdt1.m addpath(src, problems); rng(42); % 固定随机种子保证结果可复现 params struct(... popSize, 100, ... % 种群规模 gen, 200, ... % 迭代代数 dim, 30, ... % 决策变量维度ZDT1 常用 30 pc, 0.9, ... % 交叉概率 pm, 1/30, ... % 变异概率 1 / dim etaC, 20, ... % SBX 分布指数 etaM, 20, ... % 多项式变异分布指数 lb, zeros(1,30), ...% 变量下界 ub, ones(1,30)); % 变量上界 [pop, obj] nsga2(zdt1, params); plot(obj(:,1), obj(:,2), k., MarkerSize, 8); xlabel(f_1); ylabel(f_2); title(NSGA-II on ZDT1); grid on; axis tight;在 MATLAB 命令行里执行run_zdt1跑完你会看到一张散点图100 个黑色圆点从左上到右下排成一条弧线。正常实现下这些点应该形成一条“左高右低”的凸曲线左端大约在 (0, 1)右端在 (1, 0)。rng(42)固定随机流之后每次跑出的图略有不同但整体形状一致——如果两次跑出来前沿完全一样反而可疑说明种群多样性有问题。这里有一个新手容易走的弯路直接拿真实工程问题跑跑完看结果“好像还行”就收工了。正确流程应该是先在理论前沿已知的测试函数上验证确认实现没写错再换真实问题。ZDT1 的作用就是这块“试金石”。4.2 结果判读用理论前沿检验实现是否写错ZDT1 的理论前沿是一条精确曲线当 g1 时f2 1 - sqrt(f1)f1 从 0 到 1 变化。这条线是判断算法实现正确性的标尺x linspace(0, 1, 200); trueFront 1 - sqrt(x); hold on; plot(x, trueFront, r-, LineWidth, 1.5); legend(NSGA-II, ZDT1 理论前沿);把理论前沿叠上去之后你会看到三种可能的结果第一黑点贴着红线分布略有一点波动——实现正确算法正常工作第二黑点明显偏离红线成群结队地漂在红线一侧——实现有问题多半是非支配排序或拥挤度的判定逻辑写错了第三黑点只覆盖红线的中间一段两端空白——种群多样性不足需要加大 popSize 或调整变异强度。对双目标问题肉眼对比就是最强的验证手段不需要任何统计指标。如果你改代码后前沿形状变了先别急着调参回去查排序函数是不是对目标矩阵的维度处理错了。4.3 记录收敛指标用 IGD 或 HV 量化对比纯肉眼判断在写论文或做算法对比时不够需要量化指标。两个通用的选择IGD反世代距离衡量算法前沿到理论前沿的平均距离越小越好HV超体积衡量前沿覆盖的面积越大越好。常见做法是每隔若干代记录一次前沿算一次 IGD画一条收敛曲线。IGD 的 MATLAB 实现很直白function igd computeIGD(approxFront, trueFront) % approxFront: 算法得到的 Pareto 点集 % trueFront: 理论前沿采样点集 d 0; for i 1:size(trueFront, 1) dists sqrt(sum((approxFront - trueFront(i,:)).^2, 2)); d d min(dists); end igd d / size(trueFront, 1); endIGD 每一代算一次就能看出算法在 50 代、100 代、200 代各处于什么水平。实际操作中把nsga2主循环里每一代的obj输出来存到 cell 数组里跑完后统一计算比改主循环更稳妥。下面是固定 popSize100用默认参数跑 ZDT1 的典型观察结果迭代代数前沿形态IGD 大致量级50点少且零散两端明显缺失0.05 ~ 0.1200点覆盖完整贴住理论前沿0.005 ~ 0.02500与理论前沿基本重合肉眼难辨差异 0.005注意“大致量级”是测试函数上常见的参考区间不代表固定结论——IGD 受采样点数和随机种子影响。重要的是趋势代数增加IGD 应单调下降如果 200 代比 50 代反而更大说明种群在“退化”这基本可以确定代码某个环节丢了精英解。5. 避坑与调参从结果异常到性能瓶颈的排查路线5.1 结果怎么判读先看前沿形状再谈收敛指标拿到一张 Pareto 图我习惯按三步走第一步看前沿形状对不对是不是符合问题的理论特征第二步看端点有没有覆盖到两端稀疏说明边界解丢了第三步看 IGD 或 HV 的收敛曲线确认代数加大会继续改善而不是震荡。很多人在第一步就直接翻车——前沿形状不对还对着指标调参等于在错误的基础上做微调越调越偏。记住指标只能告诉你“离理论多远”不能告诉你“方向对不对”。形状错了指标再小也是自欺欺人。5.2 踩坑记录这一节写我实际踩过的几个坑每条都是“现象 → 原因 → 解决”希望能让你少走几趟夜路。坑一种群规模一上 500单代耗时肉眼可见地变慢。现象popSize 从 100 提到 500每次迭代等了快一分钟跑 200 代等于半小时起步。原因fastNonDominatedSort是双重 for 循环O(MN²) 的复杂度在 N500 时就到了瓶颈。解决把内层判定向量化——用all(pop(i,:) pop, 2)一次性算 i 与所有个体的支配关系把双重循环降为一重。向量化后 500 个体的排序时间能压到原来的十分之一。坑二拥挤度算出来一半是 NaN。现象crowdingDistance返回的数组里出现 NaN环境选择排序时这些个体被排到最后种群多样性崩了。原因某个目标维度上所有个体取值完全相同fmax - fmin 0没短路就出现除零。解决在分母计算前加if fmax fmin, continue; end。出现这一般说明目标函数写错了返回了常数。坑三跑完画图所有点全都挤在 (1,0) 附近。现象Pareto 前沿缩成一个点簇几乎看不到曲线形状。原因变异概率pm设成了 1 而不是 1/dim等于每个个体每代每个变量都在变异搜索变成了随机扰动算法退化成随机游走。解决把pm改为1/dimZDT1 的 30 维问题就是 1/30。坑四换了目标函数后rank 全变成 1。现象非支配排序结果里所有个体 rank 1选择机制失效。原因目标函数返回了列向量整个种群的目标矩阵变成了 M×N 而不是 N×Mall(pop(i,:) pop(j,:))逐元素比较时维度不匹配MATLAB 自动广播后逻辑全错。解决检查目标函数返回值一定是 1×M 的行向量且evaluateAll堆叠后严格是 N×M。坑五相同参数跑两次结果差得离谱。现象两次跑出的前沿形状大方向一致但具体点位分布毫无规律IGD 波动超过 30%。原因没设随机种子MATLAB 每次启动的随机流不同。解决主脚本开头加rng(42)或者把随机种子作为参数传进来方便批量实验时固定变量。要跑多次取统计平均时再故意换种子。5.3 调参优先级先动哪个参数才不会白跑一晚上调参最忌讳一上来就动 popSize 和 gen。我感觉的优先级是先确认前沿形状正确形状都不对直接回去改代码然后调 etaC 和 etaM这两个参数决定探索与开发的平衡——前沿覆盖不足就降 etaM 增强变异前沿震荡太大就升 etaC 让子代更稳定最后才考虑 popSize 和 gen它们是“堆资源”的参数只有在算法本身已经正常时才值得加。很多同学一上来把 popSize 加到 1000结果跑一晚上发现前沿还是缺一块因为问题的根在变异算子不在种群规模。当然如果画完图发现端点总差一截边界解那时候加 popSize 是最有效的——边界解需要足够的采样密度才能被保留下来。6. 二次开发从换测试函数到加约束的三种改造手法6.1 接口约定是整个改装的起点这套代码能改到什么程度完全由problems/目录下的函数接口决定。约定很简单输入一维行向量x输出一维行向量f。遵守这个约定任何测试问题都能无缝替换。你自己写的目标函数也得遵守第一行接收行向量最后一行返回行向量中间随便怎么写。看清楚这个接口后改造的套路就清晰了——不是在主循环里塞业务逻辑而是新写一个problems函数主循环一行都不用动。比如把 ZDT1 换成 ZDT2非凸前沿只需要新建zdt2.m把目标函数里的f2改成1 - (f1/g)^2即可nsga2(zdt2, params)直接就能跑。前沿从凸弧变成凹弧你会立刻看到黑点分布形状跟着变。6.2 从双目标换到三目标双目标看二维散点三目标看三维散点。改造点有两个目标函数输出三列绘图换成scatter3f [f1, f2, f3]; % 三目标测试函数如 DTLZ2 scatter3(obj(:,1), obj(:,2), obj(:,3), 20, filled); xlabel(f_1); ylabel(f_2); zlabel(f_3); view(135, 30); % 调整视角避免平面遮挡 grid on;核心算法部分不需要改——非支配排序、拥挤度、环境选择本来就是按任意 M 目标写的。注意拥挤度计算里for m 1:M天然支持三维唯一要检查的是你的目标函数三列的量纲是否接近量纲差太大会让拥挤度被大目标主导前沿会偏向某个目标。处理办法是对目标值做归一化后再传入排序函数或者调大那个目标的权重。6.3 加入不等式约束一个 penalty 就够用工程问题几乎都带约束。NSGA-II 原论文用的是约束支配比较入门级的简洁做法是 penalty 函数把约束违反量加权加到目标值上违反约束的解自然会被淘汰function [f, cv] constrainedZDT1(x) % 带一个不等式约束的 ZDT1 变体 f zdt1(x); c 0.5 - f(1); % 示例约束要求 f1 0.5 cv max(0, c); % 约束违反量满足约束时为 0 f f 100 * cv; % 违反约束的目标值被放大 endcv是约束违反量100是惩罚系数。常见做法是惩罚系数取 100 或 1000太小则约束形同虚设太大会让可行解区域边缘的目标值失真。运行后画图时你只会看到 f1 小于等于 0.5 的那半边前沿——约束在效果上切掉了不可行区域。如果要做得更精细再考虑约束支配比较那样需要额外改fastNonDominatedSort给每个个体加一个CV字段支配判定时先比可行性再比目标值。从那以后我每次拿到任何多目标优化代码第一件事都是先跑一遍理论前沿已知的 ZDT1确认前沿形状对上了再往上面接工程问题。这个习惯帮我挡掉了至少三回“跑出漂亮图但解全错”的翻车现场——那些图漂亮是因为代码忠实地把某个错误的逻辑跑到了收敛而不是因为问题本身该长那样。做优化的同学得先学会相信眼睛看到的形状再去信数值指标。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助?