1 数学模型
用于模拟电价的模型是一个简化形式的混合模型，如下图1所示。其根本驱动因素是天然气价格和气温。该模型在内部捕获了驱动因素与电价的关系之间的关系，以及与一天中的时间、一周中的哪一天和节假日的关系。天然气价格模型是一种简化的均值回归随机微分方程模型。温度模型是参数函数和时间序列模型的组合。
模拟的天然气和电价路径被送入调度算法，该算法计算工厂的最优调度，以产生一系列现金流，可用于计算风险现金流。
图1
部分代码：
%% 每小时温度建模和模拟 % 这个例子演示了每小时干球温度拟合一个非线性温度模型。 % 温度序列被建模为两个分量的总和，一个确定性的非线性函数， % 用来解释给定年份某一时刻的季节或预期温度，另一个随机分量，用来解释实际温度与平均值的偏差。 %% 导入数据 %数据集从先前创建的MAT-file中加载。它由一系列日期的向量和相应的历史记录温度组成。 clear load Data\TempSeries whos %% 季节性成分（确定性） % 确定性或期望的温度分量用正弦和模型建模，激励的是数据中观察到的温度和周期性的物理本质。 % 利用曲线拟合工具箱对模型参数进行估计。统计工具箱函数NLINFIT也可以用来估计模型参数。 % 计算平均值(年平均值)并将其从系列中剔除 m = mean(drybulb) drybulb0 = drybulb - m; % Fit double-sine model model = fit(dates, drybulb0, 'sin2') %% 可视化模型精度 % 用自定义图对拟合结果进行可视化分析。图中的两个轴是相连的，这极大地使数据的视觉探索。 % 还要注意，X轴上的日期tick相对于缩放和平移是动态的。 pred = model(dates) + m; res = drybulb - pred; fitPlot(dates, [drybulb pred], res); disp(['平均绝对误差: ' num2str(mean(abs(res))) ' 华氏度']); %% 分析残差中的序列相关性 % 上述图中明显的特征之一是残差不具有序列相关性。预计这是因为高于平均气温很可能跟随高于平均气温。 % 这种序列相关性可以被明确地度量。这里用AUTOCORR和PARCORR函数来显示序列的自相关和部分自相关。 figure; subplot(2,1,1); autocorr(res,50); title('随机序列的序列相关性'); subplot(2,1,2); parcorr(res(1:1000),50); %% 利用季节AR模型对随机成分进行建模 % 可以选择对随机分量建模一个均值回复漂移SDE。 % 但是，由于季节性，我们将使用带有季节性滞后的自回归模型。 % 这里可以使用MATLAB反斜杠算子代替REGRESS函数，但这不会返回置信区间。 lags = [1 2 3 4 23 24 25 47 48 49]; Xres = lagmatrix(res, lags); [beta, betaci, res2] = regress(res, Xres); disp('滞后系数和置信区间'); disp([lags' beta betaci]) %% 对回归残差进行序列相关分析 %来自回归的残差现在应该大多是序列不相关的。 figure; subplot(2,1,1); plot(dates, res2); datetick title('回归残差及其序列相关性'); subplot(2,1,2); autocorr(res2(lags(end)+1:end),50); %% 对残差拟合一个分布 % 由于残差大多是不相关的，因此可以将它们建模为独立的从适当的分布中提取的。 % 一个t位置尺度分布可以显示出很好的拟合性。 PD = fitdist(res2, 'tlocationscale'); %% 模型总结 % 温度模型现在可以定义为， % % *平均气温' m ' % *正弦模型'模型' % *回归参数'β ' % *自相关滞后'滞后' % *剩馀概率分布PD % *可选样本前数据( 回归温度的最后一次观测值 ) tempModel = struct('m', m, 'sinmodel', model, 'reglags', lags, 'regbeta', beta, 'dist', PD, 'presample', res(end-lags(end)+1:end)); save SavedModels\TemperatureModel.mat -struct tempModel clearvars -except tempModel dates drybulb %% 仿真模型 % 我们现在可以对这个模型进行2009年的模拟，并将模拟值与2009年的观测数据进行比较。 % 注意我们已经捕获了我们模型中的相关特征。 newDates = dates(end-365*24+1:end); simTemp = simulateTemperature(tempModel, newDates, 1); %% 可视化仿真结果 ax1 = subplot(2,1,1); plot(newDates, drybulb(end-365*24+1:end)) title('实际温度'); ax2 = subplot(2,1,2); plot(newDates, simTemp); title('模拟温度'); linkaxes([ax1 ax2], 'x'); dynamicDateTicks([ax1 ax2], 'linked');
2 运行结果
应用程序中包含三个 Excel 电子表格。 PlantRiskDeployed.xlsm 是与已部署应用程序一起使用的接口。它包括调用使用 MATLAB Builder Ex 创建的 COM DLL 的已部署方法的 VBA 代码。部署项目 PlantRisk.prj 包括构建已部署应用程序所需的相关文件。这可以在 MATLAB 中打开并构建以创建 Excel 应用程序所需的 DLL。 vbamodule.bas 中的 VBA 代码包含已成为 PlantRiskDeployed 一部分的 VBA 代码，用于从 Excel 自动调用 DLL 函数。 PlantRiskSpLink.xlsm 是与 Spreadsheet Link EX 一起使用的接口，用于在 MATLAB 的运行会话中调用 MATLAB 函数。如果您不需要与不是 MATLAB 用户的其他人共享应用程序，或者在创建已部署的应用程序之前对界面进行原型设计，则此版本非常有用。 PlantRiskTemplate.xlsx 是没有集成到 MATLAB 的简单接口。 vbamodule.bas 中的 VBA 代码以及 MATLAB Builder EX 生成的 VBA 代码可以导入此文件，使其成为与 PlantRiskDeployed.xlsm 相同的功能齐全的部署应用程序.
使用界面：单击 Run Backtest 按钮调用 backtestPlantPortfolio 函数（在 MATLAB 中或在 MATLAB 生成的 DLL 中）。此函数加载工作表上指定日期范围内的历史电力和天然气价格，并计算每个发电机的每日最佳调度，从而在尊重最小运行时间限制的同时最大化每日发电利润。回测的结果是每个工厂的一组运营统计数据，例如已实现利润、运营天数、工厂运行时间百分比以及每个运营日的平均小时数。该函数还计算投资组合的总预期利润和日期范围内每个日期的现金流（利润）图。回测结果的屏幕截图如下所示。单击 Run Simulation 按钮会调用 simulationPlantPortfolio 函数（在 MATLAB 中或在 MATLAB 生成的 DLL 中）。该函数加载预训练的天然气价格、温度和电价模型，并联合模拟指定日期范围内天然气和电价的多条路径。可以在工作表上指定模拟次数。然后针对每个模拟价格路径为每个工厂执行最佳调度，以计算每个工厂以及整个投资组合的利润和运营统计数据。然后将这些平均并显示在工作表上。此外，还计算了每个资产和投资组合的风险现金流量度。该函数还创建了现金流（收益）分布的直方图，其中突出显示了风险指标
3 Matlab代码实现
部分代码：
%% 每小时温度建模和模拟 % 这个例子演示了每小时干球温度拟合一个非线性温度模型。 % 温度序列被建模为两个分量的总和，一个确定性的非线性函数， % 用来解释给定年份某一时刻的季节或预期温度，另一个随机分量，用来解释实际温度与平均值的偏差。 %% 导入数据 %数据集从先前创建的MAT-file中加载。它由一系列日期的向量和相应的历史记录温度组成。 clear load Data\TempSeries whos %% 季节性成分（确定性） % 确定性或期望的温度分量用正弦和模型建模，激励的是数据中观察到的温度和周期性的物理本质。 % 利用曲线拟合工具箱对模型参数进行估计。统计工具箱函数NLINFIT也可以用来估计模型参数。 % 计算平均值(年平均值)并将其从系列中剔除 m = mean(drybulb) drybulb0 = drybulb - m; % Fit double-sine model model = fit(dates, drybulb0, 'sin2') %% 可视化模型精度 % 用自定义图对拟合结果进行可视化分析。图中的两个轴是相连的，这极大地使数据的视觉探索。 % 还要注意，X轴上的日期tick相对于缩放和平移是动态的。 pred = model(dates) + m; res = drybulb - pred; fitPlot(dates, [drybulb pred], res); disp(['平均绝对误差: ' num2str(mean(abs(res))) ' 华氏度']); %% 分析残差中的序列相关性 % 上述图中明显的特征之一是残差不具有序列相关性。预计这是因为高于平均气温很可能跟随高于平均气温。 % 这种序列相关性可以被明确地度量。这里用AUTOCORR和PARCORR函数来显示序列的自相关和部分自相关。 figure; subplot(2,1,1); autocorr(res,50); title('随机序列的序列相关性'); subplot(2,1,2); parcorr(res(1:1000),50); %% 利用季节AR模型对随机成分进行建模 % 可以选择对随机分量建模一个均值回复漂移SDE。 % 但是，由于季节性，我们将使用带有季节性滞后的自回归模型。 % 这里可以使用MATLAB反斜杠算子代替REGRESS函数，但这不会返回置信区间。 lags = [1 2 3 4 23 24 25 47 48 49]; Xres = lagmatrix(res, lags); [beta, betaci, res2] = regress(res, Xres); disp('滞后系数和置信区间'); disp([lags' beta betaci]) %% 对回归残差进行序列相关分析 %来自回归的残差现在应该大多是序列不相关的。 figure; subplot(2,1,1); plot(dates, res2); datetick title('回归残差及其序列相关性'); subplot(2,1,2); autocorr(res2(lags(end)+1:end),50); %% 对残差拟合一个分布 % 由于残差大多是不相关的，因此可以将它们建模为独立的从适当的分布中提取的。 % 一个t位置尺度分布可以显示出很好的拟合性。 PD = fitdist(res2, 'tlocationscale'); %% 模型总结 % 温度模型现在可以定义为， % % *平均气温' m ' % *正弦模型'模型' % *回归参数'β ' % *自相关滞后'滞后' % *剩馀概率分布PD % *可选样本前数据( 回归温度的最后一次观测值 ) tempModel = struct('m', m, 'sinmodel', model, 'reglags', lags, 'regbeta', beta, 'dist', PD, 'presample', res(end-lags(end)+1:end)); save SavedModels\TemperatureModel.mat -struct tempModel clearvars -except tempModel dates drybulb %% 仿真模型 % 我们现在可以对这个模型进行2009年的模拟，并将模拟值与2009年的观测数据进行比较。 % 注意我们已经捕获了我们模型中的相关特征。 newDates = dates(end-365*24+1:end); simTemp = simulateTemperature(tempModel, newDates, 1); %% 可视化仿真结果 ax1 = subplot(2,1,1); plot(newDates, drybulb(end-365*24+1:end)) title('实际温度'); ax2 = subplot(2,1,2); plot(newDates, simTemp); title('模拟温度'); linkaxes([ax1 ax2], 'x'); dynamicDateTicks([ax1 ax2], 'linked'); %% 每小时温度建模和模拟 % 这个例子演示了每小时干球温度拟合一个非线性温度模型。 % 温度序列被建模为两个分量的总和，一个确定性的非线性函数， % 用来解释给定年份某一时刻的季节或预期温度，另一个随机分量，用来解释实际温度与平均值的偏差。 %% 导入数据 %数据集从先前创建的MAT-file中加载。它由一系列日期的向量和相应的历史记录温度组成。 clear load Data\TempSeries whos %% 季节性成分（确定性） % 确定性或期望的温度分量用正弦和模型建模，激励的是数据中观察到的温度和周期性的物理本质。 % 利用曲线拟合工具箱对模型参数进行估计。统计工具箱函数NLINFIT也可以用来估计模型参数。 % 计算平均值(年平均值)并将其从系列中剔除 m = mean(drybulb) drybulb0 = drybulb - m; % Fit double-sine model model = fit(dates, drybulb0, 'sin2') %% 可视化模型精度 % 用自定义图对拟合结果进行可视化分析。图中的两个轴是相连的，这极大地使数据的视觉探索。 % 还要注意，X轴上的日期tick相对于缩放和平移是动态的。 pred = model(dates) + m; res = drybulb - pred; fitPlot(dates, [drybulb pred], res); disp(['平均绝对误差: ' num2str(mean(abs(res))) ' 华氏度']); %% 分析残差中的序列相关性 % 上述图中明显的特征之一是残差不具有序列相关性。预计这是因为高于平均气温很可能跟随高于平均气温。 % 这种序列相关性可以被明确地度量。这里用AUTOCORR和PARCORR函数来显示序列的自相关和部分自相关。 figure; subplot(2,1,1); autocorr(res,50); title('随机序列的序列相关性'); subplot(2,1,2); parcorr(res(1:1000),50); %% 利用季节AR模型对随机成分进行建模 % 可以选择对随机分量建模一个均值回复漂移SDE。 % 但是，由于季节性，我们将使用带有季节性滞后的自回归模型。 % 这里可以使用MATLAB反斜杠算子代替REGRESS函数，但这不会返回置信区间。 lags = [1 2 3 4 23 24 25 47 48 49]; Xres = lagmatrix(res, lags); [beta, betaci, res2] = regress(res, Xres); disp('滞后系数和置信区间'); disp([lags' beta betaci]) %% 对回归残差进行序列相关分析 %来自回归的残差现在应该大多是序列不相关的。 figure; subplot(2,1,1); plot(dates, res2); datetick title('回归残差及其序列相关性'); subplot(2,1,2); autocorr(res2(lags(end)+1:end),50); %% 对残差拟合一个分布 % 由于残差大多是不相关的，因此可以将它们建模为独立的从适当的分布中提取的。 % 一个t位置尺度分布可以显示出很好的拟合性。 PD = fitdist(res2, 'tlocationscale'); %% 模型总结 % 温度模型现在可以定义为， % % *平均气温' m ' % *正弦模型'模型' % *回归参数'β ' % *自相关滞后'滞后' % *剩馀概率分布PD % *可选样本前数据( 回归温度的最后一次观测值 ) tempModel = struct('m', m, 'sinmodel', model, 'reglags', lags, 'regbeta', beta, 'dist', PD, 'presample', res(end-lags(end)+1:end)); save SavedModels\TemperatureModel.mat -struct tempModel clearvars -except tempModel dates drybulb %% 仿真模型 % 我们现在可以对这个模型进行2009年的模拟，并将模拟值与2009年的观测数据进行比较。 % 注意我们已经捕获了我们模型中的相关特征。 newDates = dates(end-365*24+1:end); simTemp = simulateTemperature(tempModel, newDates, 1); %% 可视化仿真结果 ax1 = subplot(2,1,1); plot(newDates, drybulb(end-365*24+1:end)) title('实际温度'); ax2 = subplot(2,1,2); plot(newDates, simTemp); title('模拟温度'); linkaxes([ax1 ax2], 'x'); dynamicDateTicks([ax1 ax2], 'linked'); %% 每小时温度建模和模拟 % 这个例子演示了每小时干球温度拟合一个非线性温度模型。 % 温度序列被建模为两个分量的总和，一个确定性的非线性函数， % 用来解释给定年份某一时刻的季节或预期温度，另一个随机分量，用来解释实际温度与平均值的偏差。 %% 导入数据 %数据集从先前创建的MAT-file中加载。它由一系列日期的向量和相应的历史记录温度组成。 clear load Data\TempSeries whos %% 季节性成分（确定性） % 确定性或期望的温度分量用正弦和模型建模，激励的是数据中观察到的温度和周期性的物理本质。 % 利用曲线拟合工具箱对模型参数进行估计。统计工具箱函数NLINFIT也可以用来估计模型参数。 % 计算平均值(年平均值)并将其从系列中剔除 m = mean(drybulb) drybulb0 = drybulb - m; % Fit double-sine model model = fit(dates, drybulb0, 'sin2') %% 可视化模型精度 % 用自定义图对拟合结果进行可视化分析。图中的两个轴是相连的，这极大地使数据的视觉探索。 % 还要注意，X轴上的日期tick相对于缩放和平移是动态的。 pred = model(dates) + m; res = drybulb - pred; fitPlot(dates, [drybulb pred], res); disp(['平均绝对误差: ' num2str(mean(abs(res))) ' 华氏度']); %% 分析残差中的序列相关性 % 上述图中明显的特征之一是残差不具有序列相关性。预计这是因为高于平均气温很可能跟随高于平均气温。 % 这种序列相关性可以被明确地度量。这里用AUTOCORR和PARCORR函数来显示序列的自相关和部分自相关。 figure; subplot(2,1,1); autocorr(res,50); title('随机序列的序列相关性'); subplot(2,1,2); parcorr(res(1:1000),50); %% 利用季节AR模型对随机成分进行建模 % 可以选择对随机分量建模一个均值回复漂移SDE。 % 但是，由于季节性，我们将使用带有季节性滞后的自回归模型。 % 这里可以使用MATLAB反斜杠算子代替REGRESS函数，但这不会返回置信区间。 lags = [1 2 3 4 23 24 25 47 48 49]; Xres = lagmatrix(res, lags); [beta, betaci, res2] = regress(res, Xres); disp('滞后系数和置信区间'); disp([lags' beta betaci]) %% 对回归残差进行序列相关分析 %来自回归的残差现在应该大多是序列不相关的。 figure; subplot(2,1,1); plot(dates, res2); datetick title('回归残差及其序列相关性'); subplot(2,1,2); autocorr(res2(lags(end)+1:end),50); %% 对残差拟合一个分布 % 由于残差大多是不相关的，因此可以将它们建模为独立的从适当的分布中提取的。 % 一个t位置尺度分布可以显示出很好的拟合性。 PD = fitdist(res2, 'tlocationscale'); %% 模型总结 % 温度模型现在可以定义为， % % *平均气温' m ' % *正弦模型'模型' % *回归参数'β ' % *自相关滞后'滞后' % *剩馀概率分布PD % *可选样本前数据( 回归温度的最后一次观测值 ) tempModel = struct('m', m, 'sinmodel', model, 'reglags', lags, 'regbeta', beta, 'dist', PD, 'presample', res(end-lags(end)+1:end)); save SavedModels\TemperatureModel.mat -struct tempModel clearvars -except tempModel dates drybulb %% 仿真模型 % 我们现在可以对这个模型进行2009年的模拟，并将模拟值与2009年的观测数据进行比较。 % 注意我们已经捕获了我们模型中的相关特征。 newDates = dates(end-365*24+1:end); simTemp = simulateTemperature(tempModel, newDates, 1); %% 可视化仿真结果 ax1 = subplot(2,1,1); plot(newDates, drybulb(end-365*24+1:end)) title('实际温度'); ax2 = subplot(2,1,2); plot(newDates, simTemp); title('模拟温度'); linkaxes([ax1 ax2], 'x'); dynamicDateTicks([ax1 ax2], 'linked');
