%% 铁路监控数据共享激励机制 - 完整实验脚本
% 基于三方演化博弈和Beta信誉机制的仿真实验
% 
% 使用说明：
% 1. 确保所有函数文件在当前路径或MATLAB路径中
% 2. 直接运行此脚本即可生成所有实验结果
% 3. 所有图形会自动保存到results文件夹

clear all; close all; clc;
 params = struct();
    
    % 第一组参数（条件1）
    params(1).k1 = 1;
    params(1).C2 = 5;
    params(1).D = 80;
    params(1).C3 = 100;
    params(1).p1 = 200;
    params(1).C = 40;
    params(1).m = 50;
    params(1).a = 0.5;
    params(1).n = 50;
    params(1).k2 = 1;
    
    % 第二组参数（条件2）
    params(2).k1 = 1;
    params(2).C2 = 2;
    params(2).D = 80;
    params(2).C3 = 200;
    params(2).p1 = 200;
    params(2).C = 40;
    params(2).m = 50;
    params(2).a = 0.5;
    params(2).n = 50;
    params(2).k2 = 1;
    
    % 第三组参数（条件3）
    params(3).k1 = 1;
    params(3).C2 = 10;
    params(3).D = 100;
    params(3).C3 = 80;
    params(3).p1 = 200;
    params(3).C = 40;
    params(3).m = 50;
    params(3).a = 0.5;
    params(3).n = 100;
    params(3).k2 = 1;
    
    % 第四组参数（条件4）
    params(4).k1 = 1;
    params(4).C2 = 10;
    params(4).D = 30;
    params(4).C3 = 100;
    params(4).p1 = 200;
    params(4).C = 40;
    params(4).m = 50;
    params(4).a = 0.5;
    params(4).n = 80;
    params(4).k2 = 1;

%% 初始化
fprintf('========================================\n');
fprintf('铁路监控数据共享激励机制仿真实验\n');
fprintf('基于信誉值和智能合约的三方演化博弈\n');
fprintf('========================================\n\n');

% 创建结果保存文件夹
if ~exist('results', 'dir')
    mkdir('results');
end

% 添加函数路径
addpath(pwd);

%% 1. Beta信誉机制分析（图2）
fprintf('[1/6] 正在进行Beta信誉机制分析...\n');
beta_reputation_analysis();
saveas(gcf, 'results/figure2_beta_distribution.png');
fprintf('    图2已保存\n\n');

%% 2. 主要演化博弈仿真（图5-6）
fprintf('[2/6] 正在进行主要演化博弈仿真...\n');
run_main_simulation();
saveas(gcf, 'results/figure5_6_evolution_3d.png');
fprintf('    图5-6已保存\n\n');

%% 3. 稳定性分析
fprintf('[3/6] 正在进行稳定性分析...\n');
stability_analysis(params(1),'条件一');
stability_analysis(params(2),'条件二');
stability_analysis(params(3),'条件三');
stability_analysis(params(4),'条件四');
saveas(gcf, 'results/stability_analysis.png');
fprintf('    稳定性分析结果已保存\n\n');

%% 4. 信誉参数敏感性分析（图7-8）
fprintf('[4/6] 正在进行信誉参数k1的敏感性分析...\n');
analyze_reputation_parameter();
saveas(gcf, 'results/figure7_8_reputation_analysis.png');
fprintf('    图7-8已保存\n\n');

%% 5. 非诚实收益参数分析（图9）
fprintf('[5/6] 正在进行非诚实收益N的敏感性分析...\n');
analyze_dishonest_gain_parameter();
saveas(gcf, 'results/figure9_dishonest_gain.png');
fprintf('    图9已保存\n\n');

%% 6. 监管收益参数分析（图10）
fprintf('[6/6] 正在进行监管收益D的敏感性分析...\n');
analyze_regulation_benefit_parameter();
saveas(gcf, 'results/figure10_regulation_benefit.png');
fprintf('    图10已保存\n\n');

%% 完成
fprintf('========================================\n');
fprintf('所有实验已完成！\n');
fprintf('结果已保存在 results 文件夹中\n');
fprintf('========================================\n');

%% 辅助函数：运行主仿真
function run_main_simulation()
    % 参数设置
    params = struct();
    
    % 第一组参数（条件1）
    params(1).k1 = 1;
    params(1).C2 = 5;
    params(1).D = 80;
    params(1).C3 = 100;
    params(1).p1 = 200;
    params(1).C = 40;
    params(1).m = 50;
    params(1).a = 0.5;
    params(1).n = 50;
    params(1).k2 = 1;
    
    % 第二组参数（条件2）
    params(2).k1 = 1;
    params(2).C2 = 2;
    params(2).D = 80;
    params(2).C3 = 200;
    params(2).p1 = 200;
    params(2).C = 40;
    params(2).m = 50;
    params(2).a = 0.5;
    params(2).n = 50;
    params(2).k2 = 1;
    
    % 第三组参数（条件3）
    params(3).k1 = 1;
    params(3).C2 = 10;
    params(3).D = 100;
    params(3).C3 = 80;
    params(3).p1 = 200;
    params(3).C = 40;
    params(3).m = 50;
    params(3).a = 0.5;
    params(3).n = 100;
    params(3).k2 = 1;
    
    % 第四组参数（条件4）
    params(4).k1 = 1;
    params(4).C2 = 10;
    params(4).D = 30;
    params(4).C3 = 100;
    params(4).p1 = 200;
    params(4).C = 40;
    params(4).m = 50;
    params(4).a = 0.5;
    params(4).n = 80;
    params(4).k2 = 1;
    
    % 运行仿真
    figure('Position', [50, 50, 1500, 1000]);
    subplot_labels = {'(a)', '(b)', '(c)', '(d)'};
    for condition = 1:4
        current_params = params(condition);
        tspan = [0, 50];
        num_simulations = 20;
        
        subplot(2, 2, condition);
        hold on;
        
        for sim = 1:num_simulations
            initial_state = [rand(); rand(); rand()];
            [t, trajectory] = ode45(@(t, state) evolutionary_dynamics(t, state, current_params), ...
                                    tspan, initial_state);
            
            plot3(trajectory(:, 1), trajectory(:, 2), trajectory(:, 3), ...
                  'LineWidth', 1.5, 'Color', [0.3, 0.5, 0.8, 0.5]);
        end
        text(-0.15, 1.05, subplot_labels{condition}, ...  % 左上角位置
        'Units', 'normalized', ...                  % 使用归一化坐标
        'FontSize', 16, ...                         % 较大的字号
        'FontWeight', 'bold', ...                   % 加粗
        'VerticalAlignment', 'top', ...             % 顶部对齐
        'HorizontalAlignment', 'left');             % 左对齐
        % 设置图形属性
        xlabel('x (Regulatory department)','FontSize', 10);
        ylabel('y (Sharing Party A)','FontSize', 10);
        zlabel('z (Sharing Party B)','FontSize', 10);
        title(sprintf('Conditions %d: Evolutionary trajectory', condition));
        grid on;
        box on;
        view(45, 30);
        ax = gca; ax.GridAlpha = 0.3;
        xlim([0, 1]);
        ylim([0, 1]);
        zlim([0, 1]);
        
        % 标记平衡点
        if condition == 1
            plot3(0, 0, 0, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
            text(0, 0, 0, ' E1(0,0,0)', 'FontSize', 10);
        elseif condition == 2
            plot3(0, 1, 1, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
            text(0, 1, 1, ' E4(0,1,1)', 'FontSize', 10);
        elseif condition == 3
            plot3(1, 0, 0, 'bo', 'MarkerSize', 10, 'MarkerFaceColor', 'b');
            text(1, 0, 0, ' E5(1,0,0)', 'FontSize', 10);
        else
            plot3(1, 1, 1, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
            text(1, 1, 1, ' E8(1,1,1)', 'FontSize', 10);
        end
        
        hold off;
    end
    
    sgtitle('Simulation results of three-party evolutionary game（Four conditions）', 'FontSize', 14, 'FontWeight', 'bold');
end

function dstate_dt = evolutionary_dynamics(~, state, params)
    x = state(1);
    y = state(2);
    z = state(3);
    
    % 提取参数
    k1 = params.k1;
    k2 = params.k2;
    C2 = params.C2;
    C3 = params.C3;
    D = params.D;
    p1 = params.p1;
    C = params.C;
    m = params.m;
    a = params.a;
    n = params.n;
    
    % 复制动态方程
    dx_dt = x * (1 - x) * ((C2 * D * y * z) + D - D * y * z - C3 + C2 * z * y);
    dy_dt = y * (1 - y) * (k1 * p1 * z + C * z + m * x + a * n - n - C);
    dz_dt = z * (1 - z) * (a * n + C * y + m * x + k2 * p1 * y - n - C);
    
    dstate_dt = [dx_dt; dy_dt; dz_dt];
end