北太案例-利用北太天元竞赛版本求解2024年数学建模国赛B题第2，3问

% 定义成品合格率函数function Hgl = New_cpl(Lpj1_cpl,Lpj2_cpl,Cp_cpl,x1,x2)Hgl = ((x1 + (1 - x1)*(1 - Lpj1_cpl)) * (x2 + (1 - x2)*(1 - Lpj2_cpl)) * (1 - Cp_cpl));end

% 定义产生新决策变量的函数function [jcbl_ljjc1,jcbl_bcpjc1,jcbl_cpjc1,jcbl_cpcj1,jcbl_bcpcj1] = New_jcbl(jcbl_ljjc,jcbl_bcpjc,jcbl_cpjc,jcbl_cpcj,jcbl_bcpcj)    n = [jcbl_ljjc jcbl_bcpjc jcbl_cpjc jcbl_cpcj,jcbl_bcpcj];    r = randi([1,length(n)]);    if n(r) == 0        n(r) = 1;    else        n(r) = 0;    end    jcbl_ljjc1 = n(1,1:8);    jcbl_bcpjc1 = n(1,9:11);    jcbl_cpjc1 = n(1,12);    jcbl_cpcj1 = n(1,13);    jcbl_bcpcj1 = n(1,14:16);end

% 清空工作区clear;clc;
% 参数设置Lpj1_cpl = [0.1 0.2 0.1 0.2 0.1 0.05]; % 零配件1次品率Lpj1_jccb = [2 2 2 1 8 2]; % 零配件1的检测成本Lpj2_cpl = [0.1 0.2 0.1 0.2 0.2 0.05]; % 零配件2次品率Lpj2_jccb = [3 3 3 1 1 3]; % 零配件2的检测成本Cp_cpl = [0.1 0.2 0.1 0.2 0.1 0.05]; % 成品次品率Cp_jccb = [3 3 3 2 2 3]; % 成品检测成本B_cp_dh = [6 6 30 30 10 10]; % 用户调换损失B_cp_cj = [5 5 5 5 5 40]; % 成品拆解成本
Lpj1_gmdj = 4; % 零配件1的采购成本Lpj2_gmdj = 18; % 零配件2的采购成本Cp_zpcb = 6; % 成品装配成本Cp_scsj = 56; % 成品市场售价
% 决策参数 1表示进行该决策 0表示不进行x1 = [0 1]; % 对零配件1是否进行检测x2 = [0 1]; % 对零配件2是否进行检测x3 = [0 1]; % 是否对成品进行检测x4 = [0 1]; % 是否对不合格品进行拆解


n = 1000;  % 生产的产品数量
Ccy_cb = cell(1,6); % 用于储存成本Ccy_mincb = zeros(1,6); % 用于储存最小成本Ccy_bl = cell(1,6); % 用于储存决策变量Ccy_minbl = zeros(4,6); % 用于储存最小决策变量Sr = n * Cp_scsj; % 销售收益
% 成品的装配成本Cp_zpcb_z = n * Cp_zpcb;for i = [1 2 3 4 5 6] % 6种情况    Ccy_cb_1 = zeros(1,16); % 用于储存成本    Ccy_bl_1 = zeros(4,16); % 用于储存决策变量    num = 0;    for k1 = x1        for k2 = x2            for k3 = x3                for k4 = x4                    % 计算零件1、2的采购数量                    Lpj1_cgsl = ceil(n/(1- k1*Lpj1_cpl(i)));                    Lpj2_cgsl = ceil(n/(1- k2*Lpj2_cpl(i)));
                    % 计算零件1、2的采购成本                    Lpj1_cgcb = Lpj1_cgsl * Lpj1_gmdj;                    Lpj2_cgcb = Lpj2_cgsl * Lpj2_gmdj;
                    % 计算零件1、2的检测成本                    Lpj1_jccb_z = k1 * Lpj1_cgsl * Lpj1_jccb(i);                    Lpj2_jccb_z = k2 * Lpj2_cgsl * Lpj2_jccb(i);
                    % 成品的检测成本                    Cp_jccb_z = k3 * Cp_jccb(i) * n;
                    % 合格率                    Hgl = New_cpl(Lpj1_cpl(i),Lpj2_cpl(i),Cp_cpl(i),k1,k2);
                    % 调换损失                    Dh_sc = (n - Hgl * n) * B_cp_dh(i) * (1 - k3);
                    % 两个零配件都合格的情况下成品不合格的概率                    P_f = 1 - (1 - Lpj1_cpl(i)) * (1 - Lpj2_cpl(i)) * (1 - Cp_cpl(i));
                    %拆解后零配件的再利用概率                    P_lpj1_zjy = Lpj1_cpl(i) * (1 - P_f);                    P_lpj2_zjy = Lpj2_cpl(i) * (1 - P_f);                    P_lpj_zjy = (P_lpj1_zjy + P_lpj2_zjy) / 2; % 平均再利用率
                    % 拆解成本                    Cj_cb = k4 * (1 - Hgl) * n * B_cp_cj(i);
                    % 再利用收益，扣除的成本                    Zly_cy = k4 * (1 - Hgl) * n * P_lpj_zjy * (Lpj1_gmdj + Lpj2_gmdj);
                    % 总拆解成本                    Cj_cb = Cj_cb - Zly_cy;
                    % 总成本                    Num_cb = Cp_zpcb_z + Lpj1_cgcb + Lpj2_cgcb + Lpj1_jccb_z + Lpj2_jccb_z + Cp_jccb_z + Dh_sc + Cj_cb;
                    num = num + 1;                    Ccy_cb_1(num) = Num_cb;                    xx = [k1;k2;k3;k4];                    Ccy_bl_1(:,num) = xx;                end            end        end    end    Ccy_cb{i} = Ccy_cb_1;    Ccy_mincb(i) = min(Ccy_cb_1);    Ccy_bl{i} = Ccy_bl_1;    Ccy_minbl(:,i) = Ccy_bl_1(:,find(Ccy_cb_1 == min(Ccy_cb_1),1));end
x = 1:1:16;y1 = Ccy_cb{1};y2 = Ccy_cb{2};y3 = Ccy_cb{3};y4 = Ccy_cb{4};y5 = Ccy_cb{5};y6 = Ccy_cb{6};figureplot(x,y1,'-o',x,y2,'-o',x,y3,'-o',x,y4,'-o',x,y5,'-o',x,y6,'-o')legend('情景1','情景2','情景3','情景4','情景5','情景6')title('不同情景下各方案总成本对比')xlabel('方案编号')ylabel('总成本')
figure(2);bar(Ccy_mincb);hold on;xlabel('情景');ylabel('最小总成本');title('每种情景下的最小成本');

% 清空工作区clear;clc;
% 参数设置Cp_scsj = 200; % 成品市场售价N = 1000; % 半成品数量Lj_cpl = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]; % 零件1至8的次品率Lj_jmdj = [2, 8, 12, 2, 8, 12, 8, 12]; % 零件1至8的购买单价Lj_jccb = [1, 1, 2, 1, 1, 2, 1, 2]; % 零件1至8的检测成本
Bcp_cpl = 0.1; % 半成品次品率Bcp_zpcb = 8; % 半成品装配成本Bcp_jccb = 4;% 半成品的检测成本Bcp_cjf = 10; % 半成品拆解费用
Cp_cpl = 0.1; % 成品次品率Cp_zpcb = 8; % 成品装配成本Cp_jccb = 6; % 成品的检测成本Cp_cjf = 10; % 成品拆解费用
Cp_scsj = 200; % 成品市场售价B_cp_dh = 40; % 不合格产品调换损失
% 随机生成初始决策变量Sjs_01 = randi([0,1],1,16);%决策变量jcbl_ljjc = Sjs_01(1,1:8); % 是否检测零件1至8jcbl_bcpjc = Sjs_01(1,9:11);% 是否检测半成品1、2、3jcbl_cpjc = Sjs_01(1,12); % 是否检测成品jcbl_cpcj = Sjs_01(1,13); % 是否拆解不合格成品jcbl_bcpcj = Sjs_01(1,14:16);% 是否拆解不合格半成品1、2、3
% 定义模拟退火算法的的一些参数T0 = 10; % 初始温度 1T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是 T0maxgen = 30; % 外层迭代次数Lk = 10; % 每个温度下的迭代次数alpfa = 0.95; % 温度衰减系数% 调用函数计算初始成本Z_zb0 = Sum_cb(jcbl_ljjc,jcbl_bcpjc,jcbl_cpjc,jcbl_cpcj,jcbl_bcpcj);Min_Z_zb1 = Z_zb0; % 储存最小成本Best_jcbl = [jcbl_ljjc jcbl_bcpjc jcbl_cpjc jcbl_cpcj jcbl_bcpcj]; %储存对应决策变量RESULT_Min_Z_zb1 = zeros(maxgen,1); % 用于储存每次外循环的最小成本Min_Z_zb1for iter = 1 : maxgen % 外 循 环    for i = 1 : Lk % 内 循 环           [jcbl_ljjc1,jcbl_bcpjc1,jcbl_cpjc1,jcbl_cpcj1,jcbl_bcpcj1] = New_jcbl(jcbl_ljjc,jcbl_bcpjc,jcbl_cpjc,jcbl_cpcj,jcbl_bcpcj);           Z_zb1 = Sum_cb(jcbl_ljjc1,jcbl_bcpjc1,jcbl_cpjc1,jcbl_cpcj1,jcbl_bcpcj1);            % 新 解 更 优            if Z_zb1 < Z_zb0                jcbl_ljjc = jcbl_ljjc1;                jcbl_bcpjc = jcbl_bcpjc1;                jcbl_cpjc  = jcbl_cpjc1;                jcbl_cpcj = jcbl_cpcj1;                jcbl_bcpcj = jcbl_bcpcj1;                Z_zb0 = Z_zb1;            % 原 始 解 更 优，遵循一定概率接受            else                p = exp(-(Z_zb1 - Z_zb0)/T);                if rand(1) < p                    jcbl_ljjc = jcbl_ljjc1;                    jcbl_bcpjc = jcbl_bcpjc1;                    jcbl_cpjc  = jcbl_cpjc1;                    jcbl_cpcj = jcbl_cpcj1;                    jcbl_bcpcj = jcbl_bcpcj1;                    Z_zb0 = Z_zb1;                end            end            % 判断是否更新最优解            if Z_zb0 < Min_Z_zb1                Min_Z_zb1 = Z_zb0;                Best_jcbl = [jcbl_ljjc,jcbl_bcpjc,jcbl_cpjc,jcbl_cpcj,jcbl_bcpcj];            end    end    RESULT_Min_Z_zb1(iter) = Min_Z_zb1;     T = alpfa*T;end
figure(2);plot(RESULT_Min_Z_zb1);hold on;xlabel('外循环迭代次数');ylabel('最小总成本');title('优化目标值随迭代次数变化表');

% 定义求解成本的函数function ZCB = Sum_cb(jcbl_ljjc,jcbl_bcpjc,jcbl_cpjc,jcbl_cpcj,jcbl_bcpcj)
N = 1000; % 半成品数量Lj_cpl = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]; % 零件1至8的次品率Lj_jmdj = [2, 8, 12, 2, 8, 12, 8, 12]; % 零件1至8的购买单价Lj_jccb = [1, 1, 2, 1, 1, 2, 1, 2]; % 零件1至8的检测成本
Bcp_zpcb = 8; % 半成品装配成本Bcp_jccb = 4;% 半成品的检测成本Bcp_cjf = 10; % 半成品拆解费用
Cp_cpl = 0.1; % 成品次品率Cp_zpcb = 8; % 成品装配成本Cp_jccb = 6; % 成品的检测成本Cp_cjf = 10; % 成品拆解费用
B_cp_dh = 40; % 不合格产品调换损失
% 零件1至8的采购和检测Lj_hgl = zeros(1, 8);Lj_cg_cb = zeros(1, 8);Lj_jy_cb = zeros(1, 8);for i = 1:8    Lj_hgl(i) = 1 - Lj_cpl(i) * (1 - jcbl_ljjc(i)); % 零件合格率    Lj_cg_sl = N / (1 - Lj_cpl(i) * jcbl_ljjc(i)); % 零件采购数量    Lj_cg_cb(i) = Lj_cg_sl * Lj_jmdj(i); % 零件采购成本    Lj_jy_cb(i) = jcbl_ljjc(i) * Lj_cg_sl * Lj_jccb(i); % 零件检测成本end
% 半成品1、2、3的合格率Bcp1_hgl = prod(Lj_hgl(1:3)); % 半成品1由零件1、2、3组成Bcp2_hgl = prod(Lj_hgl(4:6)); % 半成品2由零件4、5、6组成Bcp3_hgl = prod(Lj_hgl(7:8)); % 半成品3由零件7、8组成
% 半成品检测成本Bxp_jc_cb = jcbl_bcpjc .* N .* Bcp_jccb.* 3;
% 成品的合格率Cb_hgl = Bcp1_hgl * Bcp2_hgl * Bcp3_hgl * (1 - Cp_cpl * (1 - jcbl_cpjc));
% 成品的数量pp = [(1 - Bcp1_hgl*jcbl_bcpjc(1)),(1 - Bcp2_hgl*jcbl_bcpjc(2)),(1 - Bcp2_hgl*jcbl_bcpjc(2))];Cp_n = min(pp) * N;
% 成品检测成本Cp_jc_cb = jcbl_cpjc * Cp_n * Cp_jccb;
% 半成品装配成本Bcp_zp_cb = 3 * N * Bcp_zpcb;
% 成品装配成本Cp_zp_cb = Cp_n * Cp_zpcb;
% 不合格品的调换成本Bhgp_dg_cb = (1 - jcbl_cpcj) * Cp_n * (1 - Cb_hgl) * B_cp_dh;
% 拆解后再利用率p_f_b1 = sum(Lj_hgl(1:3))/3;p_f_b2 = sum(Lj_hgl(4:6))/3;p_f_b3 = sum(Lj_hgl(7:8))/2;
p_f = (Bcp1_hgl + Bcp2_hgl + Bcp3_hgl)/3;
% 不合格半成品的拆解成本Bhgbp1_cj_cb = jcbl_bcpcj(1) * N * (1 - Bcp1_hgl) * Bcp_cjf;Bhgbp2_cj_cb = jcbl_bcpcj(2) * N * (1 - Bcp2_hgl) * Bcp_cjf;Bhgbp3_cj_cb = jcbl_bcpcj(3) * N * (1 - Bcp3_hgl) * Bcp_cjf;
% 回收收益Bhbp_hc_cy = jcbl_bcpcj(1) * N * (1 - Bcp1_hgl) * sum(Lj_jmdj(1:3)) * p_f_b1 + ...    jcbl_bcpcj(2) * N * (1 - Bcp2_hgl) * sum(Lj_jmdj(4:6)) * p_f_b2 + ...    jcbl_bcpcj(3) * N * (1 - Bcp3_hgl) * sum(Lj_jmdj(7:8)) * p_f_b3;
% 不合格成品的拆解成本Bhgp_cj_cb = jcbl_cpcj * Cp_n * (1 - Cb_hgl) * Cp_cjf;
% 回收收益Bhp_hc_cy = jcbl_cpcj * Cp_n * (1 - Cb_hgl) * sum(Lj_jmdj) * p_f;
% 总成本计算ZCB = sum(Lj_cg_cb) + sum(Lj_jy_cb) + Bcp_zp_cb + sum(Bxp_jc_cb) + ...    Cp_zp_cb + Cp_jc_cb + Bhgp_dg_cb + Bhgp_cj_cb - Bhp_hc_cy + Bhgbp1_cj_cb + ...    Bhgbp2_cj_cb + Bhgbp3_cj_cb - Bhbp_hc_cy;end