% 定义参数 B = 40; c_h = 9; h = 1; lambda = 2; theta = 1; alpha

```matlab % 定义参数 B = 40; c_h = 9; h = 1; lambda = 2; theta = 1; alpha = 0.6; beta = 0.04; R = 15; PL = 4; PH = 8; muL = 0.4; % 生成Delta的数值网格 Delta_values = linspace(0, 20, 100); % 更大的范围 % 初始化q的数值网格 eps = 1e-6; % 定义一个非常小的数来避免q取1 q_values = linspace(0, 1 - eps, 100); % 确保q不取1 % 初始化存储结果的数组 results = []; % 对每个Delta值进行循环 for i = 1:length(Delta_values) current_Delta = Delta_values(i); muH = (B - current_Delta) / c_h; % 计算muH % 初始化当前Delta下的最大值和对应的q local_max = -Inf; local_max_q = 0; % 对每个q值进行循环 for j = 1:length(q_values) current_q = q_values(j); % 检查q是否满足条件 if current_q < (muH - lambda * current_q) / (lambda * (1 - current_q)) % 计算lambda1 lambda1 = lambda * (1 - current_q); % 计算S_sym对应的定积分 f_sym = @(x) 1 ./ ((1 + (theta / lambda1) - x) .* exp(lambda1 * x / theta)); S = integral(f_sym, 0, 1); % 根据q计算max值 lmbda2 = lambda * current_q + lambda * (1 - current_q) * alpha * exp(-beta * current_Delta); N_1 = (lambda1 * (muL + theta^2 * S * exp(lambda1 / theta) + muL * theta)) / (muL * theta * (1 + theta * S * exp(lambda1 / theta))); N_2 = (lmbda2 / muH) / (1 - (lmbda2 / muH)); current_max = lambda * R - h * N_1 - h * N_2; % 如果找到更大的max值，更新最大值和对应的q if current_max > local_max local_max = current_max; local_max_q = current_q; end end end % 将当前Delta下的最大值和对应的q添加到结果数组中 results = [results; current_Delta, local_max, local_max_q]; end % 输出结果 disp('Delta, Max Value, Optimal q'); disp(results); ```