ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

算法与嵌入式控制中的离散化:从坐标压缩到差分方程

算法与嵌入式控制中的离散化:从坐标压缩到差分方程 做算法题和搞嵌入式控制的人迟早都会撞上“离散化”这个词。信息学里它叫坐标离散化说白了就是把值域极大但数量稀少的点重新映射成一段紧凑的连续编号而控制领域里的离散化是把连续时间的微分方程、传递函数换成计算机能执行的差分方程。名字一样内核也相通——都是在资源有限的前提下把“没法直接算的东西”翻译成“能算的东西”。这篇博文就把这两个方向放在一起讲透给算法方向的读者提供可直接抄的模板也给写PID、调数字电源的朋友补上离散化建模和实现的完整链路。1. 内容整体设计与思路拆解1.1 离散化到底在解决什么问题先从最朴素的角度理解“离散化”三个字。一个班50个学生学号可能是20250001到20250050这种大跨度编号但老师平时点名只认1到50号。这种“把大编号按顺序重编成紧凑编号”的操作就是坐标离散化的本质。在算法题里常见的困境是这样的题目给了你 n 个操作每个操作涉及一个区间 [l, r]l 和 r 的数据范围可以到 10^9但你总不能开一个长度为 10^9 的数组去模拟吧内存直接爆炸。可是你再仔细看一眼真正会出现在题目里的坐标最多也就是 2n 个。10^9 的范围里只有 10^5 个有效点中间全是没用的空白。离散化就是把这些有效点提出来排序、去重、映射成 1、2、3……这样的连续编号然后把所有操作都映射到新坐标系里完成。数据范围从 10^9 压缩成 10^5 级别数组开得下二分查找、树状数组、线段树全部能用上。而在控制工程里离散化解决的是另一个问题单片机或者DSP本质上只能做加法、乘法和延时它不认识微分符号 d/dt。你拿一个连续PID公式 u(t)Kpe(t)Ki∫e(t)dtKd*de(t)/dt 去烧录进芯片芯片根本不知道怎么求积分和微分。你需要把积分换成求和把微分换成差分把传递函数从 s 域翻译到 z 域。这一步翻译就是控制系统里的离散化。两边对比一下你会发现思维模型惊人地一致都存在一个“大而连续”或“大而稀疏”的世界而我们只能操作“小而离散”的世界于是就需要一套规则来做映射。1.2 两种离散化的技术路线对比很多人以为“离散化”只有一个意思其实它至少分成两条技术路线对应不同的应用场景。我在下面列个表方便大家对号入座。分支核心操作典型输出主要应用领域坐标离散化值域压缩排序、去重、二分查找紧凑编号数组算法竞赛、数据结构、数据库索引连续系统离散化差分近似、Z变换、采样保持差分方程、离散传递函数自动控制、数字电源、电机控制、信号处理两条路线的共同点是“变换后信息不丢失或损失可控”。坐标离散化不会丢失任何有效点的相对顺序只是把空白区间折叠了而控制系统离散化则必须考虑采样周期和离散化方法带来的误差这个误差是可以分析和修正的。这也是为什么我在写这篇文章时一定要把两个方向并列放在一起讲。因为很多初学者在算法题里学会了坐标离散化转头去做数字电源看到“传递函数离散化”又懵了还有些做嵌入式的人反过来看算法里的离散化觉得那是“竞赛专用技巧”其实本质上都是同一个思想。理解了这种映射思维你以后碰到任何需要“压缩”或“数字化”的问题都会多一条思路。2. 坐标离散化的核心流程与代码模板2.1 四步流程排序、去重、编号、查询坐标离散化的实现其实非常固定我建议直接记成一个模板。整个流程就四步排序、去重、二分编号、用编号操作。第一步把题目里所有涉及到的坐标值全部收集到一个数组里包括区间端点、查询点、障碍物位置一个都不能少。第二步对这个数组做排序。第三步相邻重复元素只保留一个也就是去重。第四步对于任意一个原始坐标 x用二分查找在去重后的数组里找到它所在的位置这个位置的下标就是它离散化后的新编号。C 模板长这样#include bits/stdc.h using namespace std; using ll long long; // 将原始坐标数组 v 离散化返回去重排序后的映射表 vectorll discrete(vectorll v) { sort(v.begin(), v.end()); v.erase(unique(v.begin(), v.end()), v.end()); return v; } // 查询原始坐标 x 在映射表中的编号从 1 开始 int get_id(const vectorll mp, ll x) { return int(lower_bound(mp.begin(), mp.end(), x) - mp.begin()) 1; }用的时候先构造一个 vector 存所有坐标调用 discrete 得到映射表 mp之后每次想获取某个坐标的新编号就调 get_id。我习惯让编号从 1 开始因为后面配合树状数组、差分数组时下标 0 的位置经常要空出来从 1 开始能省掉很多边界判断。如果用 Python实现思路完全一样def discrete(arr): # 排序 去重 mp sorted(set(arr)) # 建立原始值到编号的映射 id_map {x: i 1 for i, x in enumerate(mp)} return mp, id_map这里唯一要注意的是 Python 的 dict 在极端情况下有哈希冲突导致性能退化的风险但在一般竞赛或工程场景里完全够用。如果对性能极度敏感还是用 C 的 lower_bound 更稳。2.2 为什么必须是排序加去重顺序为什么不能乱很多人背模板的时候只记代码不理解为什么这两步缺一不可。我拆开讲一下。先说排序。离散化之后我们要用二分查找做编号二分查找的前提就是数组有序。如果不排序你没法用 lower_bound 找到 x 应该落在哪个位置更没法保证多个坐标之间的相对大小关系在映射后依然正确。排序保证了顺序信息不变——原来小的坐标映射后编号也小原来大的编号也大。再说去重。去重是为了保证映射是一对一的。如果数组里有重复元素比如 [1, 5, 5, 9] 不去重那么 5 这个值会占两个编号位置后续你拿编号去开数组做标记时同一个坐标可能对应两个不同的桶逻辑就乱了。而且去重之后映射表的长度减小数组占用空间也更小。unique 函数的作用是把相邻重复元素移动到数组末尾并返回新的逻辑结尾迭代器所以要配合 erase 一起用否则 vector 的 size 没变后面遍历会出问题。我见过有初学者试图用 unordered_map 直接做“原始值→编号”的映射省掉二分查找。这个思路没问题但前提是你依然要先排序去重再统一建映射。如果你只对出现过的坐标建哈希表随便来一个没在表里的查询点就查不到会直接出错。所以哈希表只是替代了二分查找这一步排序去重依然是前置条件。还有一个容易踩的坑如果你在算法题里用的是浮点数坐标比如 double 类型的点建议谨慎处理。浮点数的相等判断存在精度问题两个理论上相等的 double 可能在排序后不连续导致 unique 去重失效。我个人的做法是能用整数就用整数实在要用浮点就先统一乘上一个倍数转成定点数或者用自定义精度比较器。2.3 编号从 0 开始还是从 1 开始这个问题看起来小但影响很大。从 0 开始的好处是 lower_bound 的返回值可以直接当下标用不需要加 1从 1 开始的好处是配合树状数组、线段树、差分数组时更顺手因为这些数据结构天然希望下标从 1 开始。我的建议是做算法题时统一从 1 开始。原因有两个一是树状数组下标 0 会导致 update 死循环必须从 1 开始才能规避二是差分数组在统计答案时如果编号从 1 开始前缀和从 1 扫到 m思路更清晰。如果你想从 0 开始代码里就把 get_id 返回值少加一个 1其他逻辑对应调整。关键是全篇保持一致不要在同一个程序里一会儿 0-based 一会儿 1-based否则写到最后 debug 到怀疑人生。3. 坐标离散化的几个典型使用场景3.1 区间操作问题中的端点压缩先给一个最经典也最容易上手的例子。有 n 次操作每次给一个区间 [l, r] 内的所有位置加 1之后有 m 次查询问某个位置 p 最终的值。其中 n 和 m 都是 10^5 级别但 l、r、p 的取值范围可以到 10^9。如果你不开离散化第一反应是用差分数组创建一个长度 10^95 的数组 a每次区间加就让 a[l]、a[r1]--最后做前缀和。但 10^9 的 int 数组光内存就要 4GB直接超限。离散化的做法是把所有 l、r1、p 都收集起来排序去重后得到映射表。然后每次区间加操作映射为 [id(l), id(r1)) 区间加查询点 p 映射为 id(p) 位置取前缀和。新数组长度最多是 2n m 的量级完全可以在内存里放开。这里有一个非常关键的细节r1 这个点必须单独加入待离散化列表。因为差分数组的核心思想是“区间结束之后差值要归位”你只在 [l, r] 上操作但差分的“-1”操作落在 r1 位置。如果你只离散化 l 和 rr1 这个位置第一次出现时可能不在映射表里lower_bound 会返回一个错误位置前缀和统计直接错。我刚开始写离散化时在这个点上栽过好几次印象极其深刻。3.2 矩形面积并与扫描线扫描线求矩形面积并是离散化线段树最经典的组合。核心思路是用一条垂直于 x 轴的扫描线从左往右扫过所有矩形用线段树维护当前扫描线被矩形覆盖的 y 方向总长度扫到矩形左边界就把对应 y 区间加 1扫到右边界就减 1。这里的 y 坐标如果不离散化线段树没办法维护连续的大范围区间。而且需要注意扫描线里的线段树节点维护的是“一段 y 区间的覆盖长度”不是“单个 y 点的覆盖次数”。这两个概念的区别决定了很多人在写这道题时容易出 bug。离散化的具体做法把所有矩形的 y1、y2 收集起来排序去重后作为线段树的叶子边界。建树时每个叶子节点代表一个 y 区间 (ys[i], ys[i1])所以线段树的叶子数量是映射表长度减 1。查询总覆盖长度时直接取根节点的 len 字段。实际写代码时我习惯把 y 区间的离散化结果存到一个数组 ys 里然后用“下标 i 代表区间 (ys[i], ys[i1])”来维护覆盖长度。这样区间更新的对象是离散化后的段编号而不是原始 y 坐标线段树的区间大小全部控制在 n 的量级。整个过程下来复杂度是 O(n log n)n 是矩形数量。3.3 稀疏网格与障碍物寻路另一个很实用的场景是大规模稀疏网格的寻路。想象一个 10^9 × 10^9 的棋盘上面只有 10^5 个障碍物你要从起点走到终点。直接开二维数组是痴人说梦但如果把障碍物所在的行和列拿出来离散化就可以把巨大的空白区域压缩成有限的格子。做法是把起点、终点、所有障碍物的 x 坐标放在一起排序去重y 坐标也同理处理。然后你会发现相邻两个离散化后的 x 坐标之间可能隔着很多原始坐标但这些空白区域的通行状态完全一样没有障碍物可以合并成一条“伪行”。压缩后的网格规模变成 O(m × m) 级别其中 m 是不同坐标数量然后在这个小网格上跑 BFS 或者 A* 就行。这里要注意坐标还原。BFS 找到终点后如果需要输出路径你必须能把离散化后的坐标映射回原始坐标。所以映射表不能只存编号还得保留原始坐标值通过编号索引回原值。这就是为什么我在模板里同时返回 mp 和 get_id 的原因mp 既用来查编号也可以当作反向索引数组。3.4 与排序、二分、贪心、KMP等算法的配合离散化从来不是一个孤立存在的算法它更像是给其他算法铺路的前置步骤。排序算法、二分查找是它的底层依赖树状数组、线段树是它的常用搭档。很多贪心问题里你需要对时间区间排序后做覆盖而如果时间戳本身很大你就需要先离散化再排序。KMP 这类字符串算法虽然一般不会直接用到坐标离散化但在涉及到“模式串值域过大、需要压缩字符集”的问题里离散化一样能派上用场。我的经验是看到题目数据范围中明晃晃写着“l ≤ 10^9”或者“坐标绝对值很大”就条件反射地问自己一句——所有可能被访问到的坐标点加起来有多少个如果远小于值域那大概率就是一个离散化思路。这也解释了为什么离散化是算法竞赛里的基础技能因为它的本质是对“有限有效信息”的提炼。4. 控制系统里的离散化从微分方程到差分方程4.1 计算机只能算差分方程做过嵌入式PID控制的朋友一定深有体会你在教科书上看到的连续PID公式写得漂漂亮亮但落到单片机代码里你只能写这样的句子integral error * dt; derivative (error - last_error) / dt。这背后的数学本质就是离散化。为什么不能直接算微分因为微分是对时间的极限操作它需要连续时间上的信息而数字控制器在每一个采样周期只能看到当前时刻和过去时刻的采样值。你能得到的只有一个个离散的点最多拿这些点的差值去近似导数。于是微分方程就被替换成了差分方程连续传递函数就变成了脉冲传递函数。这里引入一个最核心的变换关系用后向差分近似微分算子 s。数学表达式是s ≈ (1 - z^{-1}) / T其中 T 是采样周期z^{-1} 表示一拍延迟。这个式子非常直观连续域的 s 是微分算子对应离散域的“当前值减上一拍的值再除以采样周期”。在C语言里就是 (error - last_error) / dt。我之所以先说后向差分是因为它在工程里最常用、最稳定、最容易实现。后向差分把 s 域左半平面映射到 z 平面的一个圆内连续系统稳定则离散系统也稳定这对写控制代码的人来说是很大的安心保障。4.2 位置式PID的离散化实现连续PID公式是这样的u(t) Kp * e(t) Ki * ∫e(t)dt Kd * de(t)/dt对它做离散化采样周期为 T把积分换成矩形法累加把微分换成后向差分得到u(k) Kp * e(k) Ki * T * Σ_{i0}^{k} e(i) Kd * (e(k) - e(k-1)) / T这个过程就是“位置式PID 用离散化差分方程”的实际落地。代码可以直接写成下面这样typedef struct { float Kp, Ki, Kd; float T; // 采样周期单位秒 float integral; // 积分累加值 float prev_error; // 上一次误差 } PID; float pid_update(PID* pid, float setpoint, float measurement) { float error setpoint - measurement; // 积分项矩形法累加 pid-integral error * pid-T; // 微分项后向差分 float derivative (error - pid-prev_error) / pid-T; pid-prev_error error; // 计算输出 return pid-Kp * error pid-Ki * pid-integral pid-Kd * derivative; }这个代码就是位置式PID的最简实现可以直接放进单片机跑。但实际工程中还要加几个细节积分限幅避免长时间误差累积导致积分饱和微分滤波因为纯粹的差分对噪声极其敏感输出限幅防止控制器输出超出执行器物理范围。积分饱和是我见过最多的坑。电机堵转或者系统饱和时误差一直存在积分项一路疯涨等误差反向时积分项还没回落系统需要很久才能恢复动态。解决办法很简单在 integral 累加后做一个 clamp限制在 [-limit, limit] 内。微分噪声放大也一样如果传感器信号毛刺多e(k)-e(k-1) 会被噪声主导Kd 稍微大一点系统就开始抖。通常的做法是给微分项串联一个低通滤波器或者直接用不完全微分PID。4.3 传递函数的离散化以数字电源为例位置式PID只是控制离散化里的一个小案例。更普遍的场景是把被控对象的传递函数 G(s) 离散成 G(z)然后在DSP里写差分方程。数字电源特别典型因为市面上的数字电源控制芯片几乎都是靠离散化后的传递函数跑环路补偿的。我举一个最简单的例子被控对象是一阶惯性环节 G(s) 1 / (1 τs)τ 是时间常数。现在要把它离散化采用后向差分 s (1 - z^{-1}) / T。代入 G(s)G(z) 1 / (1 τ * (1 - z^{-1}) / T) T / (T τ - τ * z^{-1})把它改写成一阶差分方程形式y(k) τ/(Tτ) * y(k-1) T/(Tτ) * u(k)这意味着每来一次采样新的输出 y(k) 是上一拍输出 y(k-1) 和当前输入 u(k) 的线性组合。在C语言里实现就三行// 配置T 采样周期, tau 对象时间常数 float a tau / (T tau); float b T / (T tau); // 每个采样周期中断里调用 float update(float u) { static float y_prev 0.0f; float y a * y_prev b * u; y_prev y; return y; }这里面最需要强调的是系数推导过程。很多人拿到一个离散传递函数直接照着写代码却不清楚系数从哪来。其实你只需要把 s 的替换式代入 G(s)然后化简成关于 z^{-1} 的有理式最后从 G(z) Y(z)/U(z) 交叉相乘再反变换回差分方程即可。进阶一点如果把零阶保持器ZOH也考虑进来离散化公式变成 G(z) (1 - z^{-1}) * Z(G(s)/s)。这个就需要查Z变换表或者用仿真软件帮你算系数手算会麻烦不少但更贴近真实系统行为因为DAC输出本质上就是零阶保持每个采样周期内输出保持不变直到下一拍才更新。4.4 不同离散化方法的取舍离散化方法远不止后向差分一种工程里常见的还有前向差分、双线性变换Tustin、零阶保持器ZOH。我画个表总结一下各自的特点离散化方法近似公式优点缺点适用场景前向差分s ≈ (z-1)/T实现最简单可能不稳定连续域稳定不代表离散域稳定教学示例不推荐工程使用后向差分s ≈ (1-z^{-1})/T稳定、无振荡、易实现高频段频率失真较明显数字PID、一般控制回路双线性变换(Tustin)s ≈ 2/T * (1-z^{-1})/(1z^{-1})精度高s域与z域映射关系好频率有非线性畸变需预畸变修正滤波器设计、需要高精度场合零阶保持器(ZOH)严格采样保持模型最接近真实DAC输出需查Z变换表推导较复杂数字电源、需要精确建模的场合我实际调试时发现后向差分和ZOH是使用频率最高的两个方法。普通PID控制用后向差分足够误差完全在工程允许范围内而数字电源这类对动态响应要求高的场合我会优先尝试ZOH建模然后用仿真对比阶跃响应确认误差可接受后定系数。双线性变换精度虽然好但频率范围的映射会压扁高频段直接套用可能导致截止频率偏移需要做频率预畸变。5. 实际调试中的常见问题与排查技巧5.1 坐标离散化最容易翻车的三个位置第一个翻车点前面提过就是差分数组里的 r1 忘记加进待离散化列表。这个问题隐蔽性强因为小样例可能碰巧能过数据一多就开始随机出错。我的排查习惯是写完离散化数组后把所有的 l、r、r1、查询点 p 统一打印出来看一眼确认它们全部在映射表里。第二个翻车点是 unique 和 erase 写漏。很多人只写了 unique 忘记 erase或者恰好输入里没有重复元素于是这段代码一直没被发现错误。我建议养成固定写法sort 之后直接一行 v.erase(unique(v.begin(), v.end()), v.end());不要拆开写拆开就容易漏。第三个翻车点是二分查询时越界。lower_bound 返回的是第一个大于等于 x 的迭代器如果 x 不在映射表里返回的可能是 end() 或者指向一个比 x 大的元素这时减 mp.begin() 得到的是一个错误编号。解决方案有两个一是保证所有可能查询的坐标都提前加入了待离散化数组二是查询前先做 find 判断确认存在后再 lower_bound。5.2 PID离散化实现中的典型问题PID离散化实现最典型的三个问题积分饱和、微分噪声、采样时间不固定。积分饱和的做法我前面提到了 clamp。但在具体实现时要注意限幅位置既要在 integral 累加后限幅也要在最终输出上再限幅一次因为 Kp * e 和 Kd * derivative 也可能让输出饱和。微分噪声的问题最直接的方案是给微分项加上一阶低通滤波也就是不完全微分PID。具体做法是把纯微分项视作 D(s) Kd * s * E(s)给 D(s) 串联一个低通滤波器 1/(Tf*s1)然后离散化实现这样高频噪声会被过滤掉。采样时间不固定会导致 T 变化进而 Kp、Ki、Kd 的实际作用都会漂移。我见过有人用 main 主循环跑PID循环时间被各种耗时操作拉得忽长忽短系统控制质量一塌糊涂。正确做法是把PID放进固定中断里跑采样周期由定时器保证。5.3 数字电源/控制器离散化仿真核对离散化公式算完之后我强烈建议先仿真再上硬件。最快速的做法是用 Python 的 scipy.signal 对比连续传递函数和离散传递函数的阶跃响应。import numpy as np from scipy import signal tau 0.05 # 对象时间常数单位秒 T 0.001 # 采样周期单位秒 # 连续传递函数 G(s) 1 / (1 tau*s) sys_c signal.TransferFunction([1], [tau, 1]) # 零阶保持器离散化 sys_d_zoh signal.cont2discrete(([1], [tau, 1]), T, methodzoh) print(ZOH离散化分子分母:, sys_d_zoh) # 后向差分离散化 sys_d_bwd signal.cont2discrete(([1], [tau, 1]), T, methodbilinear) print(双线性变换离散化分子分母:, sys_d_bwd) # 对比阶跃响应 t, y_c signal.step(sys_c) t_d, y_d signal.dstep(sys_d_zoh, n100)跑一遍对比图你会直观看到采样周期 T 对离散化效果的影响。T 越小离散响应越接近连续响应T 太大离散系统可能振荡甚至失稳。如果仿真阶段就发现离散化和连续模型差异巨大那多半是 T 选得不合适或者离散化方法不适合当前对象。5.4 常见问题速查表现象可能原因解决方案离散化编号后查询结果错乱查询坐标未提前加入待离散化列表把所有可能出现的坐标先收集并去重差分统计结果随机偏差r1 未加入待离散化列表单独检查端点1位置的映射unique 后数组大小没变忘记 erase使用 sort unique erase 固定三连PID 输出持续震荡微分项噪声过大或 Kd 过大加入微分滤波降低 KdPID 响应迟钝恢复慢积分饱和对 integral 做 clamp 限幅离散传递函数阶跃响应飞了采样周期 T 过大或离散化方法不合适减小 T改用 ZOH 双线性变换对比验证这张表是我实际调板子和调算法时反复用到的排错清单遇到问题可以先对着看一遍比盲猜变量靠谱得多。6. 个人经验与后续扩展建议最后分享几个我自己的使用习惯希望能帮你少走弯路。坐标离散化这边我把模板固定成了一套每次写题直接复制然后根据题目类型微调。这样省去了每次重写排序去重的精力也降低了写错基本逻辑的概率。另外我始终保留“收集所有坐标点”这一步不到万不得已不靠临时二分去现补坐标因为临时补很容易漏。控制系统离散化这边我的体会是“仿真永远先于上板”。无论公式推导得多顺畅都先拿 Python 或 MATLAB 把连续模型、离散模型、加控制器的闭环模型全跑一遍确认稳定性再动手写嵌入式代码。调数字电源时这个习惯真的救过我很多次某些极点偏移导致的问题在仿真里一眼就能看出来在硬件上可能要烧好几个板子才能发现。后续如果你对这个话题有兴趣可以继续往这几个方向扩展模型预测控制MPC里的离散预测模型卡尔曼滤波器的离散化实现以及更底层的数值积分方法梯形法、龙格-库塔法在控制器里的应用。离散化不是一道孤立的算法题它贯穿了从算法竞赛到真实工业控制的整条链路值得反复理解和练习。
返回列表