API tools faq
paste
Haxton_Sale2

oldNC_Fission_small_matlab

Haxton_Sale2
Sep 20th, 2023 (edited)
1,848
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
MatLab 11.85 KB | None | 0 0
raw download clone embed print report
  1. %Optimization for integer number of cells
  2.  
  3. Expert = true;
  4.  
  5. oxide = true;
  6. decayHastener = false;
  7. maxHeat = 9999;
  8.  
  9. enableTBU = true;
  10. enableLE = true;
  11. enableHE = true;
  12. enableMOX = false;
  13. %   Concatenate with comma(,) not semicolon(;)
  14.  
  15. % fuelList = ["HEU-235","LEU-233","HEU-233",...
  16. %     "LEN-236","MOX-241",...
  17. %     "LEP-239","HEP-239","LEP-241","HEP-241","LEA-242","HEA-242",...
  18. %     "LECm-243","HECm-243","LECm-245","HECm-245","LECm-247","HECm-247",...
  19. %     "LEB-248","HEB-248","LECf-249","HECf-249","LECf-251","HECf-251"];
  20.  
  21. fuelList = ["placeholder"];
  22.  
  23. % ["TBU","LEU-235","HEU-235","LEU-233","HEU-233",
  24. %     "LEN-236","HEN-236","MOX-239","MOX-241",
  25. %     "LEP-239","HEP-239","LEP-241","HEP-241","LEA-242","HEA-242",
  26. %     "LECm-243","HECm-243","LECm-245","HECm-245","LECm-247","HECm-247",
  27. %     "LEB-248","HEB-248","LECf-249","HECf-249","LECf-251","HECf-251"];
  28.  
  29. elementLimit = false;
  30. % if it is true, the program will try to use multiple types of elements whenever possible;
  31. % if it is false, "unbalanced" element combinations become possible
  32.  
  33. %   RF / RTG / (number) - number makes it search for a hybrid route;
  34. goal = RTG
  35.  
  36. %   Starting resources
  37. Thorium = 1800;
  38. Uranium = 1650;
  39. cycle = 23262; %Sum of all weights in a Mineralis ritual (E2E) * ticks it takes for 1 ore
  40. multiplier = 1;
  41. maxCellN = 101; %will be scaled inverse to the heat^heatFactor of the fuel
  42. reactorN = 7;
  43.  
  44. heatFactor = 0;
  45. efficiencyMax = 600;
  46. efficiencyMin = 100;
  47. effFactor = 0;
  48.  
  49. % fuelOrder = ["TBU";"LEU-235";"HEU-235";"LEU-233";"HEU-233";
  50. %     "LEN-236";"HEN-236";"MOX-239";"MOX-241";
  51. %     "LEP-239";"HEP-239";"LEP-241";"HEP-241";"LEA-242";"HEA-242";
  52. %     "LECm-243";"HECm-243";"LECm-245";"HECm-245";"LECm-247";"HECm-247";
  53. %     "LEB-248";"HEB-248";"LECf-249";"HECf-249";"LECf-251";"HECf-251"];
  54.  
  55. %%
  56.  
  57.  
  58. TBU = zeros(27,1);
  59. TBU(1) = 1;
  60. TBU = logical(TBU);
  61.  
  62. LE = zeros(27,1);
  63. LE([2:2:6, 10:2:27],1) = 1;
  64. LE = logical(LE);
  65.  
  66. HE = zeros(27,1);
  67. HE([3:2:7, 11:2:27],1) = 1;
  68. HE = logical(HE);
  69.  
  70. MOX = zeros(27,1);
  71. MOX([8, 9],1) = 1;
  72. MOX = logical(MOX);
  73.  
  74. fuelOrder = ["TBU";"LEU-235";"HEU-235";"LEU-233";"HEU-233";
  75.     "LEN-236";"HEN-236";"MOX-239";"MOX-241";
  76.     "LEP-239";"HEP-239";"LEP-241";"HEP-241";"LEA-242";"HEA-242";
  77.     "LECm-243";"HECm-243";"LECm-245";"HECm-245";"LECm-247";"HECm-247";
  78.     "LEB-248";"HEB-248";"LECf-249";"HECf-249";"LECf-251";"HECf-251"];
  79.  
  80. isoOrder = ["Th-230";"Th-232";
  81.     "U-233";"U-235";"U-238";
  82.     "Np-236";"Np-237";
  83.     "Pu-238";"Pu-239";"Pu-241";"Pu-242";
  84.     "Am-241";"Am-242";"Am-243";
  85.     "Cm-243";"Cm-245";"Cm-246";"Cm-247";
  86.     "Bk-247";"Bk-248";"Cf-249";"Cf-250";"Cf-251";"Cf-252"];
  87.  
  88. heat=[
  89.     18;     50;     300;    60;     360;                %"TBU";"LEU-235";"HEU-235";"LEU-233";"HEU-233";
  90.     36;     216;    57.5;   97.5;                       %"LEN-236";"HEN-236";"MOX-239";"MOX-241";
  91.     40;     240;    70;     420;    94;     564;        %"LEP-239";"HEP-239";"LEP-241";"HEP-241";"LEA-242";"HEA-242";
  92.     112;    672;    68;     408;    54;     324;        %"LECm-243";"HECm-243";"LECm-245";"HECm-245";"LECm-247";"HECm-247";
  93.     52;     312;    116;    696;    120;    720];       %"LEB-248";"HEB-248";"LECf-249";"HECf-249";"LECf-251";"HECf-251"
  94.  
  95. if Expert
  96.     heat = heat * 1.2;
  97. end
  98.  
  99. if oxide == true
  100.     heat([1:7, 10:27])=heat([1:7, 10:27])*1.25;
  101. end
  102.    
  103. % eff = -heat;
  104. % eff = -log(heat);
  105. eff = min(heat) ./ (heat).^effFactor + 1;
  106. eff = eff - min(eff);
  107. eff = eff ./ max(eff) * (efficiencyMax-efficiencyMin) + efficiencyMin;
  108.  
  109. eff = eff./100;
  110.  
  111. tick = [
  112.     144000; 72000;  72000;  64000;  64000;              %"TBU";"LEU-235";"HEU-235";"LEU-233";"HEU-233";
  113.     102000; 102000; 84000;  56000;                      % "LEN-236";"HEN-236";"MOX-239";"MOX-241";
  114.     92000;  92000;  60000;  60000;  54000;  54000;      % "LEP-239";"HEP-239";"LEP-241";"HEP-241";"LEA-242";"HEA-242";
  115.     52000;  52000;  68000;  68000;  78000;  78000;      % "LECm-243";"HECm-243";"LECm-245";"HECm-245";"LECm-247";"HECm-247";
  116.     86000;  86000;  60000;  60000;  58000;  58000];     % "LEB-248";"HEB-248";"LECf-249";"HECf-249";"LECf-251";"HECf-251"
  117.  
  118. if Expert
  119.     tick = tick/2;
  120. end
  121.  
  122. CC = cycle ./ tick * multiplier; %Correction matrix
  123.    
  124. RFtNormal = [
  125.     60;     120;    480;    144;    576;                % "TBU";"LEU-235";"HEU-235";"LEU-233";"HEU-233";
  126.     90;     360;    155.4;  243.6;                      % "LEN-236";"HEN-236";"MOX-239";"MOX-241";
  127.     105;    420;    165;    660;    192;    768;    % "LEP-239";"HEP-239";"LEP-241";"HEP-241";"LEA-242";"HEA-242";
  128.     210;    840;    162;    648;    138;    552;    % "LECm-243";"HECm-243";"LECm-245";"HECm-245";"LECm-247";"HECm-247";
  129.     135;    540;    216;    864;    225;    900];    % "LEB-248";"HEB-248";"LECf-249";"HECf-249";"LECf-251";"HECf-251"
  130.    
  131. if Expert
  132.     RFtNormal = RFtNormal * 6;
  133. end
  134.  
  135. RFtNormal = RFtNormal .* CC .* eff;
  136.  
  137. RFtOxide = RFtNormal*1.4;
  138. RFtOxide(8:9) = RFtNormal(8:9);
  139.  
  140. powerNormal = tick .* RFtNormal / 1000;
  141. powerOxide = tick .* RFtOxide / 1000;
  142.  
  143. separatorMatrix = [
  144.     1   9   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  145.     0   0   0   1   9   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0];
  146.  
  147. %fuelMatrix is 27*24; CC is 27*1
  148. fuelMatrix = [
  149.     0   -81 16  8   0   8   32  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  150.     0   0   0   -9  -32 0   8   0   8   8   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  151.     0   0   0   -36 -25 0   16  0   4   0   24  0   0   0   0   0   0   0   0   0   0   0   0   0
  152.     0   0   -9  0   -72 0   0   0   4   4   32  0   0   24  0   0   0   0   0   0   0   0   0   0
  153.     0   0   -36 0   -45 32  8   0   0   0   16  0   0   8   0   0   0   0   0   0   0   0   0   0
  154.     0   0   0   0   0   -9  -68 0   0   0   32  0   8   20  0   0   0   0   0   0   0   0   0   0
  155.     0   0   0   0   16  -36 -45 8   8   0   32  0   0   0   0   0   0   0   0   0   0   0   0   0
  156.     0   0   0   0   -32 0   0   0   -9  0   12  0   0   8   4   0   0   0   0   0   0   0   0   0
  157.     0   0   0   0   -32 0   0   0   0   -1  8   0   0   0   0   0   8   0   0   0   0   0   0   0
  158.     0   0   0   0   0   0   0   0   -1  0   -48 0   0   0   4   0   28  0   0   0   0   0   0   0
  159.     0   0   0   0   0   0   0   0   -36 0   -45 8   24  0   0   8   24  0   0   0   0   0   0   0
  160.     0   0   0   0   0   0   0   0   0   -9  -68 0   4   8   0   0   48  0   0   0   0   0   0   0
  161.     0   0   0   0   0   0   0   0   0   -36 -45 8   0   0   0   8   24  24  0   0   0   0   0   0
  162.     0   0   0   0   0   0   0   0   0   0   0   0   -9  -72 8   8   40  8   0   0   0   0   0   0
  163.     0   0   0   0   0   0   0   0   0   0   0   0   -36 -45 0   16  32  8   8   0   0   0   0   0
  164.     0   0   0   0   0   0   0   0   0   0   0   0   0   0   -9  0   -40 0   16  8   8   0   0   0
  165.     0   0   0   0   0   0   0   0   0   0   0   0   0   0   -36 0   -21 0   24  8   8   0   0   0
  166.     0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   -9  -72 0   40  8   4   0   0   12
  167.     0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   -36 -45 0   48  4   4   0   8   0
  168.     0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   -72 -9  20  4   0   0   8   32
  169.     0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   -45 -36 0   8   8   0   24  24
  170.     0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   -72 -9  4   0   4   56
  171.     0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   -45 -36 0   8   8   48
  172.     0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   -9  16  8   -32
  173.     0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   -36 32  16  -29
  174.     0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   -5  -12
  175.     0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   -20 3];
  176.  
  177. fuelMatrix = diag(CC)*fuelMatrix;
  178.  
  179. decayMatrix = [
  180.     1   0   0   0   -1  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  181.     0   1   0   0   0   -1  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  182.     0   0   1   0   0   0   -1  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  183.     1   0   0   0   0   0   0   -1  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  184.     0   0   0   1   0   0   0   0   -1  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  185.     0   0   0   0   0   0   1   0   0   -1  0   0   0   0   0   0   0   0   0   0   0   0   0   0
  186.     0   0   0   0   1   0   0   0   0   0   -1  0   0   0   0   0   0   0   0   0   0   0   0   0
  187.     0   0   0   0   0   0   1   0   0   0   0   -1  0   0   0   0   0   0   0   0   0   0   0   0
  188.     1   0   0   0   0   0   0   0   0   0   0   0   -1  0   0   0   0   0   0   0   0   0   0   0
  189.     0   0   0   0   0   0   0   0   1   0   0   0   0   -1  0   0   0   0   0   0   0   0   0   0
  190.     0   0   0   0   0   0   0   0   1   0   0   0   0   0   -1  0   0   0   0   0   0   0   0   0
  191.     0   0   0   0   0   0   0   0   0   1   0   0   0   0   0   -1  0   0   0   0   0   0   0   0
  192.     0   0   0   0   0   0   0   0   0   0   1   0   0   0   0   0   -1  0   0   0   0   0   0   0
  193.     0   0   0   0   0   0   0   0   0   0   0   0   0   1   0   0   0   -1  0   0   0   0   0   0
  194.     0   0   0   0   0   0   0   0   0   0   0   0   0   1   0   0   0   0   -1  0   0   0   0   0
  195.     0   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   -1  0   0   0   0
  196.     0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1   0   0   0   0   -1  0   0   0
  197.     0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1   0   0   0   0   -1  0   0
  198.     0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1   0   0   0   0   -1  0
  199.     0   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   -1];
  200.  
  201. isotopeSeparatorRF = 50; %in kRF
  202. decayHastenerRF = 50;
  203.  
  204. A = -[separatorMatrix;
  205.     fuelMatrix;
  206.     decayMatrix]';
  207. B = zeros(1,24);
  208.  
  209. if oxide == true
  210.     f1 = [ones(1,2)*isotopeSeparatorRF, -powerOxide', ones(1,20)*decayHastenerRF];
  211. else
  212.     f1 = [ones(1,2)*isotopeSeparatorRF, -powerNormal', ones(1,20)*decayHastenerRF];
  213. end
  214. f2 = [0,0,-fuelMatrix(:,22)',0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0];
  215.  
  216. if isnumeric(goal)
  217.     f = f1 + f2*1000*goal;
  218. elseif goal == "RF"
  219.     f = f1;
  220. elseif goal == "RTG"
  221.     f = f2;
  222.     %HEB-248, LECf-249, HECf-249 produces Cf-250 for these amounts
  223. end
  224.  
  225.     C0 = heat > maxHeat;
  226.     maximumN = maxCellN ./ ((heat/heat(1)).^heatFactor);
  227.     FF = zeros(27,1);
  228.     FF = logical(FF);
  229.  
  230. N = 27;
  231. C = zeros(49,1);
  232. m = 1;
  233.  
  234. %% Building Forbidden Fuel logical matrix
  235.  
  236. FF = FF | sum(fuelList == fuelOrder, 2); % Directly from fuelOrder
  237. FF = FF | (~enableTBU & TBU);
  238. FF = FF | (~enableLE & LE);
  239. FF = FF | (~enableHE & HE);
  240. FF = FF | (~enableMOX & MOX);
  241.  
  242. %% looking for Combos
  243.  
  244. % Goal: construct array Combo, which has numbers matching fuelIndex
  245. if elementLimit
  246.     % See what categories of fuels has been completely forbidden
  247.     isoNum = [1, 4, 2, 2, 4, 2, 6, 2, 4];
  248.     isoNumC = cumsum(isoNum);
  249.     fuelRange = cell(9, 1);
  250.     fuelRange{1} = find(~FF(1));
  251.     isoNum(1) = length(fuelRange{1});
  252.     for i = 2:9
  253.         fuelRange{i} = find(~FF(isoNumC(i-1)+1:isoNumC(i))) + isoNumC(i-1);
  254.         fuelRange{i} = fuelRange{i}';
  255.         isoNum(i) = length(fuelRange{i});
  256.     end
  257.    
  258.     % "Divide" reactorN with isoNum
  259.     isoNum1 = [0, 0, 0, 0, 0, 0, 0, 0, 0];
  260.     isoNum2 = isoNum;
  261.     reactorN1 = reactorN;
  262.     while sum(isoNum2)
  263.         if sum(logical(isoNum2)) > reactorN1
  264.             break;
  265.         end
  266.         reactorN1 = reactorN1 - sum(logical(isoNum2));
  267.         isoNum1 = isoNum1 + logical(isoNum2);
  268.         isoNum2 = isoNum2 - logical(isoNum2);
  269.     end
  270.     % isoNum1 is the quotient part (locked)
  271.     % isoNum2 is the remainder (variable)
  272.    
  273.     % produce variable part
  274.     isoNum_choose = [0];
  275.     if sum(isoNum2) && reactorN1
  276.         isoNum_choose = nchoosek(find(isoNum2), reactorN1);
  277.         isoNum_rep = repmat(isoNum1, length(isoNum_choose), 1);
  278.         for i = 1:size(isoNum_choose, 1)
  279.             isoNum_rep(i, isoNum_choose(i,:)) = isoNum_rep(i, isoNum_choose(i,:)) + 1;
  280.         end
  281.     else
  282.         isoNum_rep = isoNum1;
  283.     end
  284.     Combo = [];
  285.     for i = 1:size(isoNum_rep, 1)
  286.         fuelPerm = cell(9,1);
  287.         Combo1 = [0];
  288.         for j = find(isoNum_rep(i, :))
  289.             Combo2 = nchoosek(fuelRange{j}, isoNum_rep(i, j));
  290.             Combo1 = [repelem(Combo1, size(Combo2, 1), 1), repmat(Combo2, size(Combo1, 1), 1)];
  291.         end
  292.         Combo1(:, 1) = [];
  293.         Combo = [Combo;Combo1];
  294.     end
  295. else
  296.     Combo = nchoosek(find(~FF), reactorN);
  297. end
  298.  
  299.  
  300. % fuelOrder = ["TBU";"LEU-235";"HEU-235";"LEU-233";"HEU-233"; 1-5
  301. %     "LEN-236";"HEN-236";"MOX-239";"MOX-241"; 6-9
  302. %     "LEP-239";"HEP-239";"LEP-241";"HEP-241";"LEA-242";"HEA-242"; 10-15
  303. %     "LECm-243";"HECm-243";"LECm-245";"HECm-245";"LECm-247";"HECm-247";
  304. %     "LEB-248";"HEB-248";"LECf-249";"HECf-249";"LECf-251";"HECf-251"];
  305.    
  306.  
  307. %% Loop
  308.  
  309. x = zeros(49, length(Combo));
  310. fval = zeros(1, length(Combo));
  311.  
  312. options = optimoptions('intlinprog', 'Display', 'off');
  313.  
  314. ub = ones(49,1) * 99999999;
  315. ub(3:29,1) = maximumN / multiplier;
  316. lb = zeros(1, 49);
  317.  
  318. Beq = zeros(1,49);
  319. Beq(1) = Thorium;
  320. Beq(2) = Uranium;
  321.  
  322. for k = 1:(size(Combo, 1))
  323.  
  324.     FF2 = ones(27, 1);
  325.     FF2(Combo(k, :)) = 0;
  326.     C = [ones(2,1);FF2;zeros(20,1)];
  327.  
  328.     Aeq = diag(C);
  329.     if ~mod(k, 100)
  330.         disp(string(k)+'/'+string(length(Combo)))
  331.         disp('Currently trying: ')
  332.         disp(fuelOrder(Combo(k, :))')
  333.     end
  334.        
  335.     [x2, fval2, exitflag, ~] = intlinprog(f, 3:29, A, B, Aeq, Beq, lb, ub, options);
  336.     if exitflag == 1 || exitflag == 2
  337.         x(:, k) = x2;
  338.         fval(k) = fval2;
  339.     else
  340.         fval(k) = 1000 + exitflag;
  341.     end
  342. end
  343.  
  344. fval = -fval;
  345. [~, fMAX] = max(fval);
  346. Result = x(:, fMAX);
  347. Result(3:29) = Result(3:29) .* cycle ./ tick;
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!