% 基于信誉值和智能合约的铁路监控数据共享激励机制
% 三方演化博弈仿真主程序
% Author: Railway Data Sharing Research Team
% Date: 2024

clear all; close all; clc;

%% 1. 参数设置
% 根据表3设置四组参数，对应四种条件
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;

%% 2. 运行四种条件的仿真
figure('Position', [100, 100, 1200, 800]);

for condition = 1:4
    fprintf('正在仿真条件 %d...\n', condition);
    
    % 设置当前参数
    current_params = params(condition);
    
    % 时间跨度
    tspan = [0, 50];
    
    % 运行50次仿真
    num_simulations = 50;
    
    % 创建子图
    subplot(2, 2, condition);
    hold on;
    
    % 存储所有轨迹
    all_trajectories = cell(num_simulations, 1);
    
    for sim = 1:num_simulations
        % 随机初始条件
        x0 = rand();
        y0 = rand();
        z0 = rand();
        initial_state = [x0; y0; z0];
        
        % 求解微分方程
        [t, trajectory] = ode45(@(t, state) evolutionary_dynamics(t, state, current_params), ...
                                tspan, initial_state);
        
        all_trajectories{sim} = trajectory;
        
        % 绘制轨迹
        plot3(trajectory(:, 1), trajectory(:, 2), trajectory(:, 3), ...
              'LineWidth', 1.5, 'Color', [0.3, 0.5, 0.8, 0.5]);
    end
    
    % 设置图形属性
    xlabel('x (监管部门)');
    ylabel('y (共享方A)');
    zlabel('z (共享方B)');
    title(sprintf('条件 %d: 演化轨迹', condition));
    grid on;
    box on;
    view(45, 30);
    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, 'go', 'MarkerSize', 10, 'MarkerFaceColor', 'g');
        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, 'mo', 'MarkerSize', 10, 'MarkerFaceColor', 'm');
        text(1, 1, 1, ' E8(1,1,1)', 'FontSize', 10);
    end
    
    hold off;
end

sgtitle('三方演化博弈仿真结果（四种条件）', 'FontSize', 14, 'FontWeight', 'bold');

%% 3. 生成图5和图6的二维演化过程图
generate_2D_evolution_plots(params);

%% 4. 参数敏感性分析
% 4.1 信誉参数k1的影响分析
analyze_reputation_parameter();

% 4.2 非诚实收益N的影响分析
analyze_dishonest_gain_parameter();

% 4.3 监管收益D的影响分析
analyze_regulation_benefit_parameter();

fprintf('\n仿真完成！所有图形已生成。\n');

%% 演化动力学函数
function dstate_dt = evolutionary_dynamics(~, state, params)
    x = state(1); % 监管部门选择严格监管的概率
    y = state(2); % 共享方A选择诚实共享的概率
    z = state(3); % 共享方B选择诚实共享的概率
    
    % 提取参数
    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);
    
    % 共享方A
    dy_dt = y * (1 - y) * (k1 * p1 * z + C * z + m * x + a * n - n - C);
    
    % 共享方B
    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