製造プロセスの時系列データを解析する方法｜Pythonで工程ログから品質予測へ

"""Synthetic batch-log tutorial. Numeric settings are not process recipes."""
from pathlib import Path
import argparse, json, platform
import joblib
import numpy as np
import pandas as pd
import sklearn
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Ridge
from sklearn.dummy import DummyRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error

FEATURES=['temp_mean','temp_max','hold_temp_mean','hold_temp_std',
          'ramp_rate','above_threshold_min','hold_flow_mean']

def synthetic(out):
    rng=np.random.default_rng(42); rows=[]; batches=[]
    for i in range(72):
        bid=f'B{i+1:03d}'; start=pd.Timestamp('2026-01-01 09:00',tz='UTC')+pd.Timedelta(days=i)
        target=rng.uniform(455,525); flow=rng.uniform(1.4,2.6)
        elapsed=np.arange(61)
        temp=np.where(elapsed<20,25+(target-25)*elapsed/20,
             np.where(elapsed<40,target,target-(target-60)*(elapsed-40)/20))
        temp=temp+rng.normal(0,1.0,61)
        q=70+.12*(target-490)+7*(flow-2)+rng.normal(0,.8)
        for t,v in zip(elapsed,temp):
            # One isolated gap; B008 additionally has an unacceptable long gap.
            if t==12 or (i==7 and 25<=t<=31):continue
            rows.append([bid,start+pd.Timedelta(minutes=int(t)),'temperature',float(v)])
        for t in range(0,61,2):
            rows.append([bid,start+pd.Timedelta(minutes=t),'flow',float(flow+rng.normal(0,.025))])
        batches.append([bid,start,start+pd.Timedelta(minutes=60),
                        start+pd.Timedelta(minutes=20),start+pd.Timedelta(minutes=40),
                        start+pd.Timedelta(days=3),float(q)])
    logs=pd.DataFrame(rows,columns=['batch_id','timestamp','sensor','value'])
    meta=pd.DataFrame(batches,columns=['batch_id','start','end','ramp_end','hold_end','quality_available_at','quality'])
    logs.to_csv(out/'process_logs.csv',index=False);meta.to_csv(out/'batch_quality.csv',index=False)
    return logs,meta

def prepare(logs,meta,require_quality=True):
    logs=logs.copy();meta=meta.copy()
    if not {'batch_id','timestamp','sensor','value'}<=set(logs):raise ValueError('Missing log columns')
    if not {'batch_id','start','end','ramp_end','hold_end'}<=set(meta):raise ValueError('Missing batch columns')
    if logs[['batch_id','timestamp','sensor']].isna().any().any() or meta.batch_id.isna().any():raise ValueError('Missing identity/time')
    if meta.batch_id.duplicated().any():raise ValueError('Duplicate batch ID in metadata')
    if set(logs.batch_id)!=set(meta.batch_id):raise ValueError('Log/metadata batch IDs do not match')
    def parse_time(values):
        for value in values:
            if pd.isna(value) or pd.Timestamp(value).tzinfo is None:
                raise ValueError('Timezone must be explicit: Z or +09:00')
        return pd.to_datetime(values, utc=True, errors='raise')
    logs['timestamp'] = parse_time(logs.timestamp)
    for c in ['start','end','ramp_end','hold_end']:
        meta[c] = parse_time(meta[c])
    if require_quality:
        if not {'quality','quality_available_at'} <= set(meta):
            raise ValueError('Training requires quality and quality_available_at')
        meta['quality_available_at'] = parse_time(meta.quality_available_at)
    else:
        meta['quality'] = np.nan
        meta['quality_available_at'] = meta.end
    if logs.duplicated(['batch_id','sensor','timestamp']).any():raise ValueError('Duplicate sensor timestamp')
    for c in ['value']:
        logs[c]=pd.to_numeric(logs[c],errors='raise')
        if not np.isfinite(logs[c].dropna()).all():raise ValueError('Infinite sensor value')
    meta['quality']=pd.to_numeric(meta.quality,errors='raise')
    if require_quality and not np.isfinite(meta.quality).all():raise ValueError('Missing/nonfinite quality')
    if not ((meta.start<meta.ramp_end)&(meta.ramp_end<meta.hold_end)&(meta.hold_end<meta.end)&(meta.quality_available_at>=meta.end)).all():raise ValueError('Invalid phase/availability order')
    if not set(logs.sensor)<= {'temperature','flow'}:raise ValueError('Unknown sensor; adapt schema explicitly')
    return logs.sort_values(['batch_id','sensor','timestamp']),meta

def regularize(g,grid,max_gap_s=180):
    s=g.set_index('timestamp').value.dropna().sort_index()
    if len(s)<2:raise ValueError('Insufficient observed sensor data')
    if s.index[0]>grid[0] or s.index[-1]<grid[-1]:raise ValueError('Missing boundary; no extrapolation')
    if np.max((s.index[1:]-s.index[:-1]).total_seconds())>max_gap_s:raise ValueError(f'Observed gap exceeds {max_gap_s} seconds')
    # Retain off-grid observations before time interpolation; never join batches.
    merged=s.reindex(s.index.union(grid)).sort_index().interpolate(method='time',limit_area='inside')
    out=merged.reindex(grid)
    if out.isna().any():raise ValueError('Unfilled grid value')
    return out.to_numpy()

def extract(logs,meta,out,config=None,allow_empty=False):
    config = config or {"grid_seconds":60,"max_gap_seconds":180,"threshold_c":480.0}
    rows=[];errors=[];curves=[]
    for r in meta.itertuples(index=False):
        try:
            g=logs[logs.batch_id==r.batch_id]
            if ((g.timestamp<r.start)|(g.timestamp>r.end)).any():raise ValueError('Log outside batch interval')
            grid=pd.date_range(r.start,r.end,freq=pd.Timedelta(seconds=config['grid_seconds']))
            if len(grid)<4 or grid[-1]!=r.end:raise ValueError('End must align to configured grid')
            temp=regularize(g[g.sensor=='temperature'],grid,config['max_gap_seconds'])
            flow=regularize(g[g.sensor=='flow'],grid,config['max_gap_seconds'])
            minutes=(grid-r.start).total_seconds().to_numpy()/60
            ramp_end=(r.ramp_end-r.start).total_seconds()/60
            hold_end=(r.hold_end-r.start).total_seconds()/60
            ramp=minutes<ramp_end;hold=(minutes>=ramp_end)&(minutes<hold_end)
            if ramp.sum()<2 or hold.sum()<2:raise ValueError('Insufficient phase points')
            slope=np.polyfit(minutes[ramp],temp[ramp],1)[0]
            # Left-endpoint approximation; endpoint itself carries no duration.
            duration=float(np.sum(np.diff(minutes)*(temp[:-1]>=config['threshold_c'])))
            rows.append(dict(batch_id=r.batch_id,start=r.start,end=r.end,
                temp_mean=float(np.mean(temp)),temp_max=float(np.max(temp)),
                hold_temp_mean=float(np.mean(temp[hold])),hold_temp_std=float(np.std(temp[hold],ddof=1)),
                ramp_rate=float(slope),above_threshold_min=duration,hold_flow_mean=float(np.mean(flow[hold]))))
            curves.append(pd.DataFrame({'batch_id':r.batch_id,'timestamp':grid,'elapsed_min':minutes,'temperature':temp,'flow':flow}))
        except ValueError as e:errors.append({'batch_id':r.batch_id,'reason':str(e)})
    pd.DataFrame(errors,columns=['batch_id','reason']).to_csv(out/'excluded_batches.csv',index=False)
    if not rows:
        if allow_empty:
            return pd.DataFrame(columns=['batch_id','start','end','quality',*FEATURES]),pd.DataFrame(),errors
        raise ValueError('No valid batches; inspect excluded_batches.csv')
    features=pd.DataFrame(rows)
    # Every retained batch has exactly one label row; never an inner join that hides loss.
    data=features.merge(meta[['batch_id','quality','quality_available_at']],on='batch_id',how='left',validate='one_to_one',indicator=True)
    if not (data['_merge']=='both').all():raise ValueError('Unmatched retained batch')
    data.drop(columns='_merge').to_csv(out/'batch_features.csv',index=False)
    clean=pd.concat(curves,ignore_index=True);clean.to_csv(out/'regularized_logs.csv',index=False)
    return data,clean,errors

def plots(logs,clean,data,train,test,pred,out):
    ids=clean.batch_id.drop_duplicates().iloc[:3]
    fig,ax=plt.subplots(1,2,figsize=(10,3.5))
    for bid in ids:
        d=clean[clean.batch_id==bid];ax[0].plot(d.elapsed_min,d.temperature,label=bid)
    for boundary in [20,40]:
        ax[0].axvline(boundary, color='gray', linestyle=':', linewidth=1)
    ax[0].set(xlabel='Elapsed time / min',ylabel='Temperature / C',title='Synthetic batch profiles');ax[0].legend()
    bid=ids.iloc[0];start=data.set_index('batch_id').loc[bid,'start']
    raw=logs[(logs.batch_id==bid)&(logs.sensor=='temperature')]
    d=clean[clean.batch_id==bid];ax[1].scatter((raw.timestamp-start).dt.total_seconds()/60,raw.value,s=12,label='Observed')
    ax[1].plot(d.elapsed_min,d.temperature,label='1 min grid');ax[1].set(xlim=(9,15),xlabel='Elapsed time / min',ylabel='Temperature / C',title='Short gap: synthetic interpolation');ax[1].legend()
    fig.tight_layout();fig.savefig(out/'process_profiles.png',dpi=140);plt.close(fig)
    fig,ax=plt.subplots(figsize=(8,2.5));ax.scatter(train.end,np.zeros(len(train)),label='Train: labels available');ax.scatter(test.end,np.ones(len(test)),label='Test: later batches')
    unavailable=data[(data.end<test.end.min())&(~data.batch_id.isin(train.batch_id))]
    if len(unavailable):ax.scatter(unavailable.end,np.full(len(unavailable),.5),marker='x',label='Excluded: label unavailable at cutoff')
    ax.set(yticks=[],xlabel='Batch end date / UTC',title='Chronological holdout');ax.legend(fontsize=7,loc='lower left',bbox_to_anchor=(0,1.15),ncol=3);fig.autofmt_xdate();fig.tight_layout();fig.savefig(out/'process_split.png',dpi=140);plt.close(fig)
    fig,ax=plt.subplots(1,2,figsize=(9,3.5));ax[0].scatter(test.quality,pred);v=np.r_[test.quality,pred];ax[0].plot([v.min(),v.max()],[v.min(),v.max()],'k--');ax[0].set(xlabel='Synthetic measured quality / a.u.',ylabel='Predicted quality / a.u.')
    ax[1].plot(test.end,test.quality-pred,'o-');ax[1].axhline(0,color='k');ax[1].set(xlabel='Batch end date / UTC',ylabel='Residual / a.u.');fig.autofmt_xdate();fig.tight_layout();fig.savefig(out/'process_predictions.png',dpi=140);plt.close(fig)

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('--mode', choices=['train','predict'], default='train')
    parser.add_argument('--logs')
    parser.add_argument('--quality', help='Training batch metadata including measured quality')
    parser.add_argument('--batches', help='Prediction metadata: batch_id,start,end,ramp_end,hold_end')
    parser.add_argument('--model', default='process_results/model.joblib')
    parser.add_argument('--cutoff', default='2026-02-20T00:00:00Z')
    parser.add_argument('--output', default='process_results')
    parser.add_argument('--grid-seconds', type=int, default=60)
    parser.add_argument('--max-gap-seconds', type=int, default=180)
    parser.add_argument('--threshold-c', type=float, default=480.0)
    parser.add_argument('--tolerance', type=float, default=2.0,
                        help='Educational review tolerance in quality units; choose before final evaluation')
    args = parser.parse_args()
    out = Path(args.output)
    out.mkdir(parents=True, exist_ok=True)
    config = dict(grid_seconds=args.grid_seconds,max_gap_seconds=args.max_gap_seconds,threshold_c=args.threshold_c)
    if args.grid_seconds <= 0 or args.max_gap_seconds <= 0 or args.tolerance < 0:
        parser.error('Grid/gap must be positive and tolerance nonnegative')
    if args.mode == 'predict':
        if not args.logs or not args.batches:
            parser.error('Prediction requires --logs and --batches')
        # Load only a trusted model saved by this script in the same environment.
        bundle = joblib.load(args.model)
        config = bundle['config']
        logs, meta = prepare(pd.read_csv(args.logs),pd.read_csv(args.batches),require_quality=False)
        data, clean, errors = extract(logs,meta,out,config,allow_empty=True)
        result = meta[['batch_id','end']].copy()
        valid = data[FEATURES]
        outside = ((valid < pd.Series(bundle['feature_min'])) | (valid > pd.Series(bundle['feature_max']))).any(axis=1)
        predictions = data[['batch_id']].copy()
        predictions['prediction'] = bundle['model'].predict(valid) if len(valid) else np.array([])
        predictions['status'] = np.where(outside,'review_feature_range','predicted')
        # A range check is a review signal, not a validated applicability domain.
        result = result.merge(predictions,on='batch_id',how='left',validate='one_to_one')
        result['status'] = result.status.fillna('held_log_quality')
        reason = {r['batch_id']:r['reason'] for r in errors}
        result['reason'] = result.batch_id.map(reason).fillna('')
        result.to_csv(out/'new_predictions.csv',index=False)
        print(result.to_string(index=False))
        return
    if bool(args.logs) != bool(args.quality):
        parser.error('Provide both --logs and --quality, or neither for demo')
    if args.logs:
        logs,meta = pd.read_csv(args.logs),pd.read_csv(args.quality)
        mode = 'User data; physical validity requires review'
    else:
        logs,meta = synthetic(out)
        mode = 'Synthetic educational data, not material measurements'
        # Make a label-free, later batch for end-to-end prediction practice.
        newest = meta.iloc[[-1]].copy()
        newest['batch_id'] = 'NEW001'
        for c in ['start','end','ramp_end','hold_end']:
            newest[c] = newest[c] + pd.Timedelta(days=1)
        new_logs = logs[logs.batch_id == meta.iloc[-1].batch_id].copy()
        new_logs['batch_id'] = 'NEW001'
        new_logs['timestamp'] += pd.Timedelta(days=1)
        new_logs.to_csv(out/'new_logs.csv',index=False)
        newest[['batch_id','start','end','ramp_end','hold_end']].to_csv(out/'new_batches.csv',index=False)
    logs,meta = prepare(logs,meta)
    data,clean,errors = extract(logs,meta,out,config)
    cutoff = pd.Timestamp(args.cutoff)
    if cutoff.tzinfo is None:
        raise ValueError('Cutoff timezone must be explicit')
    cutoff = cutoff.tz_convert('UTC')
    train = data[(data.end<cutoff)&(data.quality_available_at<=cutoff)].sort_values('end')
    test = data[data.end>=cutoff].sort_values('end')
    if len(train)<10 or len(test)<5:
        raise ValueError('Need >=10 train and >=5 test batches; review cutoff')
    baseline = DummyRegressor(strategy='mean').fit(train[FEATURES],train.quality)
    model = make_pipeline(StandardScaler(),Ridge(alpha=1.0)).fit(train[FEATURES],train.quality)
    pred = model.predict(test[FEATURES])
    bp = baseline.predict(test[FEATURES])
    predictions = test[['batch_id','end','quality']].copy()
    predictions['prediction'] = pred
    predictions['baseline'] = bp
    predictions['absolute_error'] = np.abs(test.quality-pred)
    predictions['review_error'] = predictions.absolute_error > args.tolerance
    predictions.sort_values('absolute_error',ascending=False).to_csv(out/'test_predictions.csv',index=False)
    pd.DataFrame({'feature':FEATURES,'coefficient_in_standardized_units':model[-1].coef_}).to_csv(out/'ridge_coefficients.csv',index=False)
    versions = dict(python=platform.python_version(),numpy=np.__version__,pandas=pd.__version__,sklearn=sklearn.__version__)
    joblib.dump(dict(model=model,features=FEATURES,config=config,cutoff=str(cutoff),versions=versions,
                     feature_min=train[FEATURES].min().to_dict(),feature_max=train[FEATURES].max().to_dict()),out/'model.joblib')
    # Earlier rolling calendar windows are development diagnostics; final test remains separate.
    reports = []
    for delta in [20,10]:
        origin = cutoff-pd.Timedelta(days=delta)
        tr = data[(data.end<origin)&(data.quality_available_at<=origin)]
        va = data[(data.end>=origin)&(data.end<origin+pd.Timedelta(days=10))]
        if len(tr)>=10 and len(va)>=5:
            m = make_pipeline(StandardScaler(),Ridge(alpha=1.0)).fit(tr[FEATURES],tr.quality)
            reports.append(dict(origin=str(origin),end=str(origin+pd.Timedelta(days=10)),train=len(tr),validation=len(va),mae=float(mean_absolute_error(va.quality,m.predict(va[FEATURES])))))
    pd.DataFrame(reports,columns=['origin','end','train','validation','mae']).to_csv(out/'period_validation.csv',index=False)
    summary = dict(mode=mode,cutoff=str(cutoff),input_batches=len(meta),retained_batches=len(data),
                   excluded_log_batches=errors,excluded_rate=len(errors)/len(meta),training_batches=len(train),test_batches=len(test),
                   unavailable_label_batches=int(((data.end<cutoff)&(data.quality_available_at>cutoff)).sum()),
                   features=FEATURES,alpha_fixed=1.0,config=config,versions=versions,
                   tolerance=args.tolerance,review_error_batches=int(predictions.review_error.sum()),
                   baseline_mae=float(mean_absolute_error(test.quality,bp)),ridge_mae=float(mean_absolute_error(test.quality,pred)),
                   ridge_rmse=float(np.sqrt(mean_squared_error(test.quality,pred))))
    (out/'summary.json').write_text(json.dumps(summary,ensure_ascii=False,indent=2))
    plots(logs,clean,data,train,test,pred,out)
    print(json.dumps(summary,ensure_ascii=False,indent=2))

if __name__ == '__main__':
    main()