guanlan-src/scripts/sim_hub/pain_yawslip.py

# -*- coding: utf-8 -*-
"""偏航滑移 —— 真实机群实测。
数据源: 50 台机组的滑移事件榜 (逐台次数 / 最大滑移角 / 逐月分布)。已脱敏: 台号重编, 不含场名与厂商。
下排三种信号为合成示意, 用于教"怎么分" —— 已机器自检: 滑移轨迹须过判据第 3、4 门。"""
import os, json, numpy as np, matplotlib
matplotlib.use("Agg"); import matplotlib.pyplot as plt
from matplotlib import font_manager as fm
for c in ["PingFang SC","Hiragino Sans GB","Songti SC","STHeiti"]:
    try:
        fm.findfont(fm.FontProperties(family=c), fallback_to_default=False)
        matplotlib.rcParams["font.sans-serif"]=[c,"DejaVu Sans"]; break
    except Exception: pass
matplotlib.rcParams["axes.unicode_minus"]=False
BG="#edf7f9"; INK="#133047"; INK2="#28536b"; MUT="#61788b"; LINE="#cbe2eb"
SPEC="#0e8898"; MEAS="#2077a8"; WARN="#d48818"; BAD="#d94b43"; OK="#23875b"
rng=np.random.default_rng(7)

SRC="/Users/yuanying/wind-analytics/outputs/guojiadian_locked/偏航滑移健康榜.json"
Y=json.load(open(SRC, encoding="utf-8"))
cnt=np.array([t["全窗滑移"] for t in Y]); mx=np.array([t["最大°"] for t in Y])
mo={}
for t in Y:
    for k,v in (t.get("permonth") or {}).items(): mo[k]=mo.get(k,0)+v
months=sorted(mo); vals=[mo[m] for m in months]
order=np.argsort(-cnt)

fig=plt.figure(figsize=(14.6,8.8), facecolor=BG)
gs=fig.add_gridspec(2,3,left=.055,right=.985,top=.775,bottom=.115,wspace=.28,hspace=.60,
                    height_ratios=[1,1.02])

# ① 逐台滑移次数 (真实)
ax=fig.add_subplot(gs[0,0]); ax.set_facecolor("#fff")
col=[BAD if c>=10 else (WARN if c>0 else "#c9dbe4") for c in cnt[order]]
ax.bar(range(len(cnt)),cnt[order],color=col,width=.86)
ax.set_xlabel(f"{len(cnt)} 台机组（按次数排序，台号已重编）",color=MUT,fontsize=10)
ax.set_ylabel("全窗滑移次数",color=MUT,fontsize=10.5)
ax.set_title("① 一半机组一次没滑，问题集中在少数台",color=INK,fontsize=13,loc="left",pad=10)
ax.annotate(f"{int(cnt.max())} 次",(0,cnt.max()),xytext=(10,-4),textcoords="offset points",
            color=BAD,fontsize=12,weight="bold")
ax.text(.52,.72,f"零滑移 {int((cnt==0).sum())} 台\n有滑移 {int((cnt>0).sum())} 台\n合计 {int(cnt.sum())} 次",
        transform=ax.transAxes,fontsize=11.5,color=INK2,va="top")

# ② 最大滑移角 (真实)
ax2=fig.add_subplot(gs[0,1]); ax2.set_facecolor("#fff")
nz=mx[mx>0]
ax2.hist(nz,bins=np.arange(0,50,4),color=MEAS,alpha=.85,edgecolor="#fff")
ax2.axvline(20,color=BAD,ls="--",lw=1.8)
ax2.text(.97,.93,f"超 20° 的有 {int((mx>20).sum())} 台",transform=ax2.transAxes,ha="right",va="top",color=BAD,fontsize=11.5,weight="bold")
ax2.set_xlabel("单次最大滑移角 (°)",color=MUT,fontsize=10.5)
ax2.set_ylabel("台数",color=MUT,fontsize=10.5)
ax2.set_title(f"② 最深一次滑了 {mx.max():.1f}°",color=INK,fontsize=13,loc="left",pad=10)

# ③ 逐月 (真实) —— 有断点
ax3=fig.add_subplot(gs[0,2]); ax3.set_facecolor("#fff")
cols=[BAD if v>=40 else MEAS for v in vals]
ax3.bar(range(len(months)),vals,color=cols,width=.72)
ax3.set_xticks(range(len(months))); ax3.set_xticklabels([m[2:] for m in months],rotation=45,ha="right",fontsize=8.6)
ax3.set_ylabel("当月滑移次数",color=MUT,fontsize=10.5)
ax3.set_title("③ 不是均匀发生的，有明显断点",color=INK,fontsize=13,loc="left",pad=10)
i0=months.index("2026-03") if "2026-03" in months else None
if i0 is not None:
    ax3.annotate("此月 0 次",(i0,1),xytext=(-6,34),textcoords="offset points",
                 arrowprops=dict(arrowstyle="->",color=INK2,lw=1.2),color=INK2,fontsize=10.5)
ax3.text(.03,.93,"随后两月合计 %d 次" % (mo.get("2026-04",0)+mo.get("2026-05",0)),
         transform=ax3.transAxes,color=BAD,fontsize=11.5,weight="bold",va="top")

# 下排: 三种信号怎么分 (合成, 教学)
t=np.linspace(0,6,720)
wd=8*np.sin(2*np.pi*t/9.5)+rng.normal(0,1.1,t.size).cumsum()*0.05
pos_s=np.zeros_like(t); cur=0.
for i in range(t.size):
    cur+=0.012*(wd[i]-cur); pos_s[i]=cur
pos_s+=rng.normal(0,.05,t.size)
pos_n=np.zeros_like(t); cmd_n=np.zeros_like(t); cur=0.
for i in range(t.size):
    if abs(wd[i]-cur)>6.0: cur=wd[i]; cmd_n[i]=1
    pos_n[i]=cur
pos_n+=rng.normal(0,.12,t.size)
pos_e=wd*0.9+rng.normal(0,.15,t.size)
for j in rng.choice(np.arange(60,t.size-60),5,replace=False): pos_e[j:]+=rng.choice([-14,-9,11,16])
# 自检: 滑移须过门3(缓慢连续) 门4(与风载同侧)
assert np.corrcoef(np.gradient(pos_s),np.gradient(wd))[0,1]>0, "滑移信号未过门4"
assert np.abs(np.diff(pos_s)).max()<0.5, "滑移信号未过门3"

for k,(name,colr,desc,pos,cmd) in enumerate([
    ("滑移",BAD,"无指令 · 缓慢连续 · 跟着风载方向",pos_s,np.zeros_like(t)),
    ("正常偏航",OK,"有指令 · 越死区才动 · 阶跃式",pos_n,cmd_n),
    ("编码器故障",WARN,"位置突跳 · 与风向无关 · 可正可负",pos_e,np.zeros_like(t))]):
    a=fig.add_subplot(gs[1,k]); a.set_facecolor("#fff")
    a.plot(t,wd,color=MEAS,lw=1.4,alpha=.6,label="风向")
    a.plot(t,pos,color=colr,lw=2.3,label="机舱位置")
    if cmd.max()>0:
        for x in t[cmd>0]: a.axvline(x,color=OK,lw=.9,alpha=.45)
        a.plot([],[],color=OK,lw=.9,label="偏航指令")
    a.set_title(f"④ {name}",color=colr,fontsize=13,loc="left",pad=26,weight="bold")
    a.text(.0,1.045,desc,transform=a.transAxes,fontsize=10.4,color=INK2,va="bottom")
    a.set_xlabel("时间 (小时)",color=MUT,fontsize=10)
    if k==0: a.set_ylabel("角度 (°)",color=MUT,fontsize=10)
    a.legend(loc="upper left",fontsize=9,frameon=False,labelcolor=INK2)

for a in fig.axes:
    a.tick_params(colors=MUT,labelsize=9.2)
    for s in a.spines.values(): s.set_color(LINE)
    a.grid(color=LINE,lw=.7,alpha=.55)

fig.text(.055,.945,"偏航滑移：制动压不住风，机舱被推着走",color=INK,fontsize=20,weight="bold")
fig.text(.055,.900,f"上排是 {len(cnt)} 台机组的实测滑移事件。一半机组一次没滑，问题集中在少数台 —— "
                   "这正是机群相对筛查能做的事。",color=INK2,fontsize=12.5)
fig.text(.055,.862,"下排是三种在信号上长得很像的现象。判滑移必须四门同时成立："
                   "无偏航指令 · 风速够大 · 变化缓慢连续 · 方向与风载同侧。四门缺一就会误判。",color=MUT,fontsize=11)
fig.text(.055,.022,"上排为实测事件统计（台号已重编，不含场名与厂商）；下排为合成示意信号，"
                   "用于教学「怎么分」，其滑移轨迹已机器自检通过判据第 3、4 门。",color=MUT,fontsize=9.3)
out=os.path.join(os.path.dirname(__file__),"assets","pain_yawslip.png")
fig.savefig(out,dpi=150,facecolor=BG); print("SAVED",out)
print(f"真实: {len(cnt)} 台 | 滑移 {int(cnt.sum())} 次 | 零滑移 {int((cnt==0).sum())} 台 | 最大角 {mx.max():.1f}° | >20° {int((mx>20).sum())} 台")