A petrochemical production scheduling optimization method based on decomposition coordination and topology analysis

Abstract

The application provides a petrochemical production scheduling optimization method based on decomposition coordination and topology analysis, and belongs to the technical field of calculation, and specifically comprises the following steps: obtaining petrochemical device scheduling instance data and constructing a discrete time petrochemical production scheduling global model, and dividing the scheduling problem corresponding to the global model into a plurality of production device sub-problems and a blending and delivery sub-problem according to a production topology structure; then performing linear relaxation initialization to obtain an initial feasible upper bound, performing decomposition coordination iteration among the sub-problems through an alternating direction multiplier method, continuously recording the integer mode trajectory of each device in the iteration process, and identifying a stable integer skeleton according to the stability; finally, constructing a restricted mixed integer linear programming repair model based on a candidate fixed skeleton, obtaining a feasible scheduling scheme satisfying all the constraints of the global model, and outputting the same. Through the scheme, the solving difficulty of a large-scale petrochemical production scheduling problem can be effectively reduced, and the generation efficiency and engineering executability of the scheduling scheme are improved.

Description

Technical Field

[0001] This invention relates to the field of computing technology, and in particular to a method for optimizing petrochemical production scheduling based on decomposition coordination and topology analysis. Background Technology

[0002] Currently, among the existing petrochemical production scheduling optimization methods, centralized mixed integer programming models have a large variable scale and long direct solution time at fine time granularity. Although methods based on Lagrange decomposition can reduce the scale of a single solution, there is still a strong coupling relationship between the material supply and demand and operation mode among the subproblems after decomposition, making coordination difficult. Moreover, the decomposition results are often difficult to directly form a feasible solution that meets the complete production constraints, and a lot of subsequent repair work is still required.

[0003] It is evident that there is an urgent need for a petrochemical production scheduling optimization method based on decomposition coordination and topology analysis that offers both high optimization efficiency and adaptability. Summary of the Invention

[0004] In view of this, embodiments of the present invention provide a petrochemical production scheduling optimization method based on decomposition coordination and topology analysis, which at least partially solves the problems existing in the prior art.

[0005] This invention provides a petrochemical production scheduling optimization method based on decomposition coordination and topology analysis, comprising: Step 1: Obtain the production units, discrete time periods, operating modes, material flow, inventory, product blending, and order delivery data of the petrochemical production system to be scheduled, form a production scheduling instance, and establish a production topology based on the connection relationships between production units, materials, virtual storage tanks, product blending links, and order delivery links. Step 2: Based on the production scheduling example, construct a discrete-time global model for petrochemical production scheduling. The global model aims to minimize the comprehensive cost within the scheduling cycle. At the production scheduling level, it describes the constraints between the operating status of production units, mode switching process, material input and output, inventory changes, product blending, and order delivery. The global model is a mixed-integer linear programming model. Step 3: Based on the production topology, the scheduling problem corresponding to the global model is divided into multiple production unit sub-problems and one blending and delivery sub-problem. Each production unit sub-problem corresponds to one production unit, and the blending and delivery sub-problem retains the component oil inventory, finished oil inventory, blending flow rate, blending ratio, quality specifications, and order delivery constraints. Step 4: Transform the cross-subsystem connection relationships in the global model into coordination residuals, wherein the coordination residuals include material-oil balance residuals, component-oil connection residuals, and coupling device mode residuals, which are used to quantify the degree of deviation of the solutions of each subproblem from the cross-boundary coupling relationships. Step 5: Perform linear relaxation initialization, relax the integer variables in the global model into continuous variables and solve the linear relaxation model to obtain a reference solution that reflects the global scheduling trend. After rounding the corresponding binary variables in the reference solution, perform the first round of feasibility repair on the integer trajectory of the fixed entry device to obtain the initial feasible upper bound. Step 6: Perform decomposition coordination iteration guided by the alternating direction multiplier method. In each round of iteration, solve the subproblems of each production unit and the coordination and delivery subproblems respectively. Calculate the original residual of this round according to the definition of coordination residual, update the coordination objective, penalty parameters and Lagrange multipliers, and perform multi-round iterative coordination. Step 7: During the iterative process of the alternating direction multiplier method, continuously record the integer pattern trajectories returned by the subproblems of each production unit, and generate candidate fixed skeletons based on the stability of the steady-state pattern variables of each unit in the most recent multiple iterations. Step 8: Construct a restricted mixed integer linear programming repair model based on the candidate fixed skeleton, fix the integer trajectories of some stable devices, and re-optimize the integer variables of the unfixed devices as well as all continuous logistics, inventory, re-mixing and delivery variables to obtain a feasible scheduling scheme that satisfies all constraints of the global model. Step 9: Repeat steps 6 to 8 until the termination condition is met, and output the final scheduling scheme.

[0006] According to a specific implementation of the present invention, the production device subproblem retains the device's own steady-state mode variables, transition start-up variables, transition occupancy variables, feed boundary, yield relationship and local production cost; The transition start-up variable is used to indicate whether the production unit has started the switching process from the source operating mode to the target operating mode, and the transition occupancy variable is used to indicate whether the production unit is in a transition occupancy state during the switching from the source operating mode to the target operating mode.

[0007] According to a specific implementation of an embodiment of the present invention, the material oil balance residual The expression is: ; in, Indicates production material oil A collection of devices, Indicates the consumption of material oil A collection of devices, Indication device During the period Input flow, Indication device During the period For material oil ; output flow rate; The component oil connection residual The expression is: ; ; in, Indicates the production of component oils A collection of devices, This indicates the time period for reconciliation and delivery. For component oil Input volume, Indication device During the period Producing component oil Traffic, and These represent the component oils in the virtual tank. During the period and time period Inventory levels Indicates the blending and delivery side component oils The change in inventory of virtual containers between adjacent time periods; The coupling device mode residual The expression is: ; in, and Representing the devices and device During the period Is it in mode? .

[0008] According to a specific implementation of the present invention, the coordination target is updated using a historical smoothing method, and the update formula is: ; in, Indicates the reconciliation target index. This indicates the updated coordination objective. Indicates the first The coordination goal of the wheel Indicates historical weight, Indicates by the first The original coordination value obtained from the solution to the wheel problem.

[0009] According to a specific implementation of an embodiment of the present invention, in step 6, the penalty parameter is adaptively updated based on the relative relationship between the original residual and the dual residual. Let the first The original residual norm and the dual residual norm are respectively and The residual ratio is defined as: ,when Located in the preset balance range Internal time maintenance ,in The imbalance threshold parameter; when When the value exceeds the preset balance range, the penalty parameter is updated as follows: ; in, To update the index proportionally, To smooth out the weights, and These are the limiting coefficients for single-wheel increase and single-wheel decrease, respectively. and These are the lower and upper bounds of the penalty parameter, respectively. Indicates will Restricted to closed intervals Inside.

[0010] According to a specific implementation of the present invention, the Lagrange multiplier is updated as follows: ; in, Indicates the first The coordination residuals that are next in line Indicates the first The Lagrange multipliers of the wheel, Indicates the first The Lagrange multipliers of the wheel have not yet undergone scalar adjustment; If the penalty parameter changes, the multiplier is adjusted for synchronization scaling: ; in, and They represent the first Round iteration and the first Wheel penalty parameters, This indicates the number obtained after scale adjustment. The Lagrange multipliers of the wheel; If the penalty parameter does not change, then .

[0011] According to a specific implementation of an embodiment of the present invention, step 7 specifically includes: Suppose the system is in the most recent The trajectory of the steady-state variable recorded in the round of iteration is For the device , its first wheel stability Defined as: ; in, Indication device The set of steady-state variables, Indicates the index of the steady-state variable. The cardinality of a set is the number of elements contained in the set. This is an indicator function that takes the value 1 if the condition within the parentheses is true, and 0 otherwise. For devices with the same stability, the device trajectory is selected first to enter the candidate fixed frame in order of increasing original residual value.

[0012] According to a specific implementation of an embodiment of the present invention, the processing procedure of the restricted mixed-integer linear programming repair model includes: If the current candidate skeleton is different from the previous skeleton, then perform a restricted mixed integer linear programming repair based on the current candidate skeleton, and compare the repaired objective function value with the historical best upper bound. If it is better than the historical best upper bound, then update the historical best feasible solution and the corresponding objective value. If the current candidate skeleton is exactly the same as the previous skeleton, skip this round of repair.

[0013] According to a specific implementation of an embodiment of the present invention, the termination condition is: The historical optimal upper bound has not improved within several consecutive rounds, the integer trajectories of each device tend to stabilize, the candidate skeleton appears repeatedly, or the improvement of the original residual enters a plateau period, or the preset maximum number of iterations is reached.

[0014] The beneficial effects of this invention are as follows: By combining production topology decomposition, ADMM subproblem coordination iteration, and stable integer skeleton repair, this invention addresses the problems of existing methods, such as limited efficiency in directly solving the complete model, complex coordination of coupling relationships after decomposition, and difficulty in converting local results into a complete feasible scheduling scheme. This method reduces the scale of a single solution by utilizing production topology decomposition, coordinates coupling relationships across subproblems such as materials, inventory, and operating modes using ADMM subproblem coordination iteration, and improves the complete feasibility of the final scheme through stable integer skeleton and restricted mixed integer programming repair. This enhances the solution efficiency and executability of the scheduling scheme for short-term scheduling problems in large-scale petrochemical production. Attached Figure Description

[0015] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0016] Figure 1A flowchart illustrating a petrochemical production scheduling optimization method based on decomposition coordination and topology analysis, provided in an embodiment of the present invention; Figure 2 A factory topology diagram provided for embodiments of the present invention; Figure 3 This is a sub-problem decomposition diagram provided in the embodiments of the present invention; Figure 4 This is an algorithm framework diagram provided for an embodiment of the present invention. Detailed Implementation

[0017] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0018] The following specific examples illustrate the implementation of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. The present invention can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that, in the absence of conflict, the following embodiments and features in the embodiments can be combined with each other. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0019] It should be noted that various aspects of embodiments within the scope of the appended claims are described below. It will be apparent that the aspects described herein can be embodied in a wide variety of forms, and any particular structure and / or function described herein is merely illustrative. Based on this invention, those skilled in the art will understand that one aspect described herein can be implemented independently of any other aspect, and two or more of these aspects can be combined in various ways. For example, any number of aspects set forth herein can be used to implement the device and / or practice the method. Additionally, this device and / or method can be implemented using structures and / or functionalities other than one or more of the aspects set forth herein.

[0020] It should also be noted that the illustrations provided in the following embodiments are only schematic representations of the basic concept of the present invention. The illustrations only show the components related to the present invention and are not drawn according to the actual number, shape and size of the components in the actual implementation. In the actual implementation, the form, quantity and proportion of each component can be arbitrarily changed, and the layout of the components may also be more complex.

[0021] Furthermore, specific details are provided in the following description to facilitate a thorough understanding of the examples. However, those skilled in the art will understand that the described aspects can be practiced without these specific details.

[0022] This invention provides a petrochemical production scheduling optimization method based on decomposition coordination and topology analysis, which can be applied to the petrochemical production scheduling optimization process in industrial production scenarios.

[0023] See Figure 1 This is a flowchart illustrating a petrochemical production scheduling optimization method based on decomposition coordination and topology analysis, provided by an embodiment of the present invention. Figure 1 and Figure 4 As shown, the method mainly includes the following steps: Step 1: Obtain the production units, discrete time periods, operating modes, material flow, inventory, product blending, and order delivery data of the petrochemical production system to be scheduled, form a production scheduling instance, and establish a production topology based on the connection relationships between production units, materials, virtual storage tanks, product blending links, and order delivery links. In practice, data on production units, discrete time periods, operating modes, material flow, inventory, product blending, and order delivery of the petrochemical production system to be scheduled are obtained to form a production scheduling instance.

[0024] Establish a production topology based on the connections between production units, materials, virtual storage tanks, product blending processes, and order delivery processes, such as... Figure 2 As shown, the input and output materials of each device, the production and consumption devices of materials, the upstream and downstream relationships between devices, and the connection relationship between the production device and the blending and delivery link are determined.

[0025] Step 2: Based on the production scheduling example, construct a discrete-time global model for petrochemical production scheduling. The global model aims to minimize the comprehensive cost within the scheduling cycle. At the production scheduling level, it describes the constraints between the operating status of production units, mode switching process, material input and output, inventory changes, product blending, and order delivery. The global model is a mixed-integer linear programming model. In practical implementation, a discrete-time petrochemical production scheduling basic model is constructed based on the petrochemical production scheduling example formed in step one. The model describes the constraints between the operating status of production units, mode switching process, material input and output, inventory changes, product blending and order delivery at the production scheduling level, without elaborating on the reaction mechanism inside each production unit.

[0026] Regarding the operation of production facilities, steady-state mode variables are set. Over-occupancy variables and transition start variables ; in, Used to indicate production equipment During the period Is it in steady-state operation mode? ; Used to indicate production equipment During the period Is it in running mode? To running mode The transitional occupancy state during switching; Used to indicate production equipment Is it during the time period? Startup from running mode To running mode The switching process. All the above variables are binary variables.

[0027] For any production unit and any time period, establish mutually exclusive constraints on operating states, ensuring that the production unit is in either a steady-state operating mode or a transitional occupancy state within the same time period. The mutually exclusive constraints on operating states are as follows: ; in, A collection of production facilities. A set of discrete time periods For production unit indexing, Index for discrete time periods and These represent the source running mode and the target running mode, respectively.

[0028] Based on the fixed duration required for switching production units between different operating modes, establish the correlation constraints between transition start-up variables, transition occupancy variables, and steady-state mode variables before and after the switch: ; ; ; ; in, Indicates production equipment By operating mode Switch to running mode The required fixed duration This represents the last discrete time period within the scheduling cycle. This indicates the time period index used for summation. The above constraints are used to define the source operating mode before mode switching, the target operating mode after mode switching, the continuous transition period corresponding to a single switching process, and to ensure that the same switching process cannot be restarted before it has ended. By setting transition start variables and transition occupancy variables respectively, the start time and continuous occupancy period of mode switching are distinguished.

[0029] Regarding material input and output, set steady-state input variables. and transitional input variables ; in, Indicates production equipment During the period In steady-state operation mode Input amount at time; Indicates production equipment During the period In running mode To running mode Input values ​​during the transition between states.

[0030] The steady-state input and transitional input satisfy the corresponding upper and lower limits of the processing capacity, and the output of the production unit to different intermediate materials is determined based on the steady-state material yield and the transitional material yield.

[0031] The output is: ; ; in, For intermediate material collection, For intermediate material indexing; Indicates production equipment During the period For intermediate materials The output quantity; Indicates production equipment In steady-state operation mode Down to intermediate materials The yield; Indicates production equipment In the operating mode To running mode During the transitional state of switching, intermediate materials The yield.

[0032] For material oil that is directly transferred between production units without passing through virtual storage tanks, a material supply and demand balance constraint is established within the same time period to match the output of the production unit that generates the corresponding material oil with the input of the production unit that consumes the material oil.

[0033] The material supply and demand balance constraint is: ; in, For material oil collection For material oil index; Indicates the production of material oil A collection of production facilities; Indicates the consumption of material oil A collection of production facilities; Indicates production equipment During the period For material oil The output quantity; Indicates production equipment During the period The amount of input.

[0034] For component oils and refined oils, inventory recursion relationships are established separately. Component oil inventory is updated based on the previous period's inventory, the current period's inbound inventory, and blending consumption; refined oil inventory is updated based on the previous period's inventory, the current period's blending input, and order delivery volume. The inventory recursion relationship is as follows: ; ; in, It is a collection of component oils. It is a collection of refined oil products. For order collection; , and These are component oils, finished oils, and order indexes, respectively. and They represent the component oils respectively. During the period and time period Inventory levels; and They represent refined oil products. During the period and time period Inventory levels; Indicates production equipment During the period To component oil The quantity of goods received; Indicates component oil During the period Outbound quantity used for product blending; Indicates refined oil During the period The amount of goods entering the warehouse as a result of the blending process; Indicates time period Inbound orders Delivery of refined oil products The quantity.

[0035] The inventory levels of component oils and finished oils meet their respective upper and lower limits. The input quantities of the production unit under steady-state operation and transitional occupancy conditions meet their respective upper and lower limits of processing capacity.

[0036] In terms of product blending, blending ratio constraints are established based on the permissible usage ratios of each component oil in the finished oil, and product quality constraints are established based on the oil property parameters of each component oil and the property specifications of the finished oil, so that the blended finished oil meets the corresponding component ratio requirements and quality specifications.

[0037] Regarding order delivery, order delivery constraints are established based on the demand volume, allowed start time, and deadline for delivery for each order.

[0038] The order delivery quantity is only allowed to take values ​​within the corresponding delivery time window, and the cumulative delivery quantity shall not exceed the demand quantity of the corresponding order; for order demands that have not been completed by the end of the scheduling cycle, a stockout penalty shall be included in the objective function.

[0039] The discrete-time petrochemical production scheduling basic model aims to minimize the overall cost within the scheduling cycle. This overall cost includes crude oil cost, transitional operating cost of production units, steady-state operating cost of production units, inventory holding cost, and order shortage penalty. The objective function is: ; in, Indicates the time period of the atmospheric distillation unit ATM crude oil imports, Indicates the price of crude oil; Indicates production equipment By operating mode To running mode Transition operating costs during switching; Indicates production equipment In steady-state operation mode Operating costs below; This represents the inventory holding cost coefficient; Indicates order The out-of-stock penalty coefficient; Indicates order For refined oil Demand; This represents the undelivered quantity of the corresponding order. The complete petrochemical production scheduling model is modeled as a MILP model.

[0040] Step 3: Based on the production topology, the scheduling problem corresponding to the global model is divided into multiple production unit sub-problems and one blending and delivery sub-problem. Each production unit sub-problem corresponds to one production unit, and the blending and delivery sub-problem retains the component oil inventory, finished oil inventory, blending flow rate, blending ratio, quality specifications, and order delivery constraints. In specific implementation, such as Figure 3 As shown in the dashed box, in order to reduce the difficulty of solving the complete model, the system divides the complete scheduling problem into multiple production unit sub-problems and a scheduling and delivery sub-problem based on the production process topology.

[0041] Each production unit subproblem corresponds to one unit, preserving the unit's own steady-state mode, transitional start-up, transitional occupation, feed boundary, yield relationship, and local production cost; the blending and delivery subproblem preserves the component oil inventory, finished oil inventory, blending flow rate, blending ratio, quality specifications, and order delivery constraints.

[0042] set up This represents production-side decision-making, including mode variables, switching variables, feed rate, and output for each production unit; Let... This represents the blending and delivery decisions, including component oil inventory, refined oil inventory, blending flow rate, and order delivery volume. The objective function of the complete problem can be summarized as: ; in, Indication device Production and operating costs and mode switching costs, This represents inventory holding costs, blending-related costs, and non-delivery penalties.

[0043] Through the above division, the system breaks down the large-scale centralized model into several smaller sub-problems, but the material connections and pattern consistency across subsystems in the original model still need to be handled through a coordination mechanism.

[0044] Step 4: Transform the cross-subsystem connection relationships in the global model into coordination residuals, wherein the coordination residuals include material-oil balance residuals, component-oil connection residuals, and coupling device mode residuals, which are used to quantify the degree of deviation of the solutions of each subproblem from the cross-boundary coupling relationships. In practice, the system transforms the cross-subsystem connection relationships in the complete model into coordinated residuals, which serve as the conflict metric and coordination basis in the ADMM iteration.

[0045] The coordination residuals include three categories: material oil balance residuals, component oil connection residuals, and coupling device mode residuals.

[0046] The first category is the residual of material-oil balance. For material-oil... and time period Since this type of material does not have virtual tanks, the output of the upstream unit on the production side should be consistent with the consumption of the downstream unit at the same time.

[0047] The system measures its deviation Defined as: ; in, Indicates production material oil A collection of devices, Indicates the consumption of material oil A collection of devices, Indication device During the period Input flow, Indication device During the period For material oil The output flow rate. When When the value approaches zero, it indicates that the supply and demand relationship of the material oil across different units is becoming more coordinated.

[0048] The second category is component oil bonding residuals. The component oils produced by the production unit enter the virtual tanks on the blending and delivery side.

[0049] After decomposition and solution, the production side provides the component oil output, while the blending and delivery side provides the component oil input demand. A buffer is allowed between the two through inventory changes. The system will handle this connection deviation. Defined as: ; in, Indicates the production of component oils A collection of devices, This indicates the time period for reconciliation and delivery. For component oil Input volume, Indication device During the period Producing component oil Traffic, This represents the change in inventory in virtual tanks for blending and delivery components between adjacent time periods, when an inventory projection mechanism is used. This represents the amount of connection adjustment absorbed by the inventory buffer within the upper and lower bounds of the inventory. This residual is used to describe the degree of coordination between production-side output, supply-side demand, and the inventory buffer.

[0050] The third category is coupled device mode residuals. For device pairs with mode coupling relationships... The system requires both to maintain pattern consistency within the same time period. This is to quantify the deviation of this consistency relationship after decomposition and solution. ,definition: ; in, and Representing the devices and device During the period Is it in mode? .

[0051] when When the value approaches zero, it indicates that the mode selection of the coupling device tends to be consistent.

[0052] The system combines the above three types of residuals into the first... Total residual vector of round iteration : ; This vector reflects the degree to which the solutions to each subproblem satisfy the cross-boundary coupling relationship, and serves as the basis for subsequent multiplier updates, penalty parameter adjustments, and stable skeleton selection.

[0053] Step 5: Perform linear relaxation initialization, relax the integer variables in the global model into continuous variables and solve the linear relaxation model to obtain a reference solution that reflects the global scheduling trend. After rounding the corresponding binary variables in the reference solution, perform the first round of feasibility repair on the integer trajectory of the fixed entry device to obtain the initial feasible upper bound. In practice, before entering the ADMM coordination iteration, the system first performs linear relaxation initialization. This initialization step preserves the continuous logistics constraints, inventory constraints, reconciliation constraints, order constraints, and schema logic relationships in the complete model, only converting integer variables... , and Relax the variables to be continuous, such that they satisfy: ; By solving the linear relaxation model, the system obtains a reference solution that reflects the global scheduling trend.

[0054] The reference solution itself may not satisfy integer feasibility, but it can provide a preliminary direction for the device's operating mode and switching structure. The corresponding binary variables in the reference solution are rounded to the nearest integer, thereby providing hot start information for subsequent restricted MILP repair.

[0055] Subsequently, the system performs the first round of feasibility repair. The system selects the inlet device as the root device and fixes its reference integer trajectory. This is because the operating structure of the crude oil processing inlet device affects the feed and component oil supply of subsequent devices, and has a strong transmission effect on the overall scheduling structure. In this embodiment, the ATM is selected as the root device.

[0056] The issues identified in the first round of repairs are as follows: ; in, Represents the complete model variables. Represents the feasible region of the complete model. Represent the complete objective function. Represents the integer trajectory of the ATM. Represents the ATM integer trajectory obtained from the solution of the subproblem, including , and The goal of this step is to obtain a feasible upper bound in a relatively short time. This information is used as a warm-up for subsequent repair processes, rather than proving global optimality in the first round.

[0057] Step 6: Perform decomposition coordination iteration guided by the alternating direction multiplier method. In each round of iteration, solve the subproblems of each production unit and the coordination and delivery subproblems respectively. Calculate the original residual of this round according to the definition of coordination residual, update the coordination objective, penalty parameters and Lagrange multipliers, and perform multi-round iterative coordination. In practice, after initialization, the system enters the ADMM-guided coordination iteration. Given the... Lagrange multipliers of the wheel Penalty parameters To achieve the coordination objective, the system solves the subproblems for each production unit, as well as the blending and delivery subproblems. To express the coordinating effect of residual information on the objective function, the system constructs the following augmented Lagrangian function. : ; in, To coordinate residual vectors across subsystems, For multiplier vectors, For penalty parameters, This indicates the transpose operation. This represents the L2 norm.

[0058] During actual execution, the system allocates coordination objectives to corresponding sub-problems, enabling parallel solution of each device's sub-problems. In each iteration, the system assigns sub-problems from production devices such as ATM, VDU, FCCU, HDS, ETH, MTBE, HTU1, HTU2, and RF, as well as blending and delivery sub-problems, to multiple threads for solution. Each production device sub-problem retains its own switching logic, feed boundaries, and yield relationships.

[0059] After all subproblems are returned, the system calculates the original residual for this round according to the residual definition in step four. Simultaneously, the system observes the changes in the coordination objective between adjacent iterations. Let the th iteration be... The target vector for wheel coordination is Then the dual residual is defined as: ; in, Indicates the first Wheel dual residuals, This represents the coordination objective vector from the previous round. The original residual reflects the degree of deviation of the current subproblem solution from the cross-boundary constraints, while the dual residual reflects the degree of change of the coordination center in adjacent rounds.

[0060] To mitigate oscillations caused by drastic changes in the coordination target, the system employs a historical smoothing approach to update the coordination target. For the [missing information]... The coordination objective is updated using the following formula: ; in, This indicates the updated coordination objective. This indicates the objectives of the previous round of coordination. Indicates historical weight, This represents the original coordination value calculated from the solution to the current wheel problem. For continuous logistics connections, It can be obtained by weighted average of the flow rates on both sides of the connection; for component oil connection, the system combines inventory projection to correct the output on the production side and the input on the blending side; for mode coupling connection, the system suppresses the repeated oscillation of mode variables between different units through historical smoothing.

[0061] The system adaptively adjusts the penalty parameter based on the relative relationship between the original residual and the dual residual. Let the first... The original residual norm and the dual residual norm are respectively and Further define the residual ratio: ; in, This represents the ratio of the original residual to the dual residual. When Located in the preset balance range Inside, keep ,in Imbalance threshold parameter. When When the value exceeds this range, the system updates the penalty parameter as follows: ; in, To update the index proportionally, To smooth out the weights, and These are the limiting coefficients for single-wheel increase and single-wheel decrease, respectively. and These are the lower and upper bounds of the penalty parameter, respectively. Indicates will Restricted to closed intervals The above update rule applies when the original residual is significantly larger than the dual residual, i.e. At that time, the update will improve This strengthens the constraint on the consistency of the subsystem; when the dual residual is significantly larger than the original residual, i.e. At that time, the update will decrease This reduces the impact of changes in the coordination objective on local subproblems; when the two types of residuals are at similar levels, the system maintains... The penalty parameter remains unchanged to avoid frequent fluctuations between adjacent iterations.

[0062] Subsequently, the system updates the multipliers: ; in, Indicates the first The coordination residuals that are next in line Indicates the first The Lagrange multipliers of the wheel, Indicates the first The Lagrange multipliers in this round have not yet undergone scaling adjustments. This update process passes the current inconsistencies across subsystems to the next iteration, thereby gradually guiding the subproblems to form a more consistent logistics and pattern structure.

[0063] If the penalty parameter changes, the multiplier is adjusted for synchronization scaling: ; If the penalty parameter does not change, then , This indicates the number obtained after scale adjustment. The Lagrange multipliers of the round. This treatment can reduce abrupt changes in multiplier terms caused by variations in the penalty parameter, keeping the multiplier terms and the quadratic penalty term consistent on a numerical scale.

[0064] It should be noted that the ADMM in this application is mainly used as a heuristic coordination mechanism. Since the scheduling model contains binary variables and fixed-duration switching logic, it does not satisfy the convexity condition on which the classic ADMM depends. Therefore, this application does not take strict convergence of ADMM to the global optimum as a necessary premise. Instead, it utilizes the residuals generated by its iterative process, the coordination objective, and the stability of integer trajectories to provide high-quality structural information for subsequent restricted MILP repair.

[0065] Step 7: During the iterative process of the alternating direction multiplier method, continuously record the integer pattern trajectories returned by the subproblems of each production unit, and generate candidate fixed skeletons based on the stability of the steady-state pattern variables of each unit in the most recent multiple iterations. In practice, during the ADMM iteration process, the system continuously records the integer pattern trajectories returned by the subproblems of each production unit. Steady-state mode variables. Transition start variables and over-occupancy variables Together they determine the operating structure of the device, among which steady-state mode variables It is the main basis for judging whether the operating framework of the device tends to be stable.

[0066] Suppose the system is in the most recent The trajectory of the steady-state variable recorded in the round of iteration is For the device , its first wheel stability Defined as: ; in, Indication device The set of steady-state variables, Indicates the index of the steady-state variable. The cardinality of a set is the number of elements contained in the set. This is an indicator function; it takes the value 1 if the condition within the parentheses is true, and 0 otherwise. Stability. The closer it is to 1, the better the device. In recent The closer the steady-state mode selection in each iteration is to a fixed state, the better. The system generates candidate fixed skeletons based on stability and residual level. For devices with high stability, the system considers their integer operating structures to have high reliability; for devices with the same stability, the system prioritizes selecting device trajectories with smaller original residuals into the candidate skeleton to avoid fixing structures that still have large cross-boundary conflicts. In this way, the system transforms the dynamic information in ADMM iterations into integer structure priors for constrained MILP repair.

[0067] Step 8: Construct a restricted mixed integer linear programming repair model based on the candidate fixed skeleton, fix the integer trajectories of some stable devices, and re-optimize the integer variables of the unfixed devices as well as all continuous logistics, inventory, re-mixing and delivery variables to obtain a feasible scheduling scheme that satisfies all constraints of the global model. In practice, the system constructs a restricted MILP repair model based on the obtained candidate fixed skeletons. Unlike directly solving the complete MILP, this repair model fixes the integer trajectories of some stable devices and only re-optimizes the integer variables of unfixed devices, as well as all continuous logistics, inventory, blending, and delivery variables. By reducing the number of free integer variables, the system can narrow the branch-bound search space and improve the efficiency of repairing the globally feasible solution.

[0068] For devices that enter the candidate skeleton, the system fixes its steady-state mode variables. Transition start variables and over-occupancy variables The reference values ​​are used; for devices that do not enter the candidate skeleton, the system retains the adjustable space of their integer variables. The repaired model still retains the relevant constraints in the complete model, so the obtained solution can meet the global scheduling feasibility requirements.

[0069] If the current candidate skeleton is different from the previous skeleton, the system performs a restricted MILP repair; if the objective function value obtained by the repair is better than the historical best upper bound UB, the historical best feasible solution and the corresponding objective value are updated.

[0070] If the current candidate skeleton is exactly the same as the one in the previous round, the system skips the duplicate repair to avoid consuming extra computation time without new structural information.

[0071] Step 9: Repeat steps 6 to 8 until the termination condition is met, and output the final scheduling scheme.

[0072] In practice, the system checks the termination conditions after each round of ADMM coordination and restricted MILP repair. When the historical optimal upper bound has not improved significantly within several consecutive rounds, the integer trajectories of each device tend to stabilize, candidate skeletons reappear, or the improvement of the original residuals enters a plateau, the system determines that the new information that can be provided by continued coordination is limited and terminates the iteration. If the maximum number of iterations is reached, the system also stops calculation, returns the current historical optimal feasible solution, and outputs the final scheduling scheme.

[0073] Following the steps outlined above, this application first transforms the petrochemical plant scheduling problem into a discrete-time mixed-integer optimization model. Then, based on the production topology, it decomposes the problem into multiple production plant subproblems and a coordination and delivery subproblem. Next, it uses ADMM to coordinate residuals, coordination objectives, multipliers, and penalty parameters to iteratively update the description of conflicts and coordination directions between subsystems. Furthermore, it identifies the stable operating skeleton and restores the decomposed coordination results to a globally feasible scheduling scheme through restricted MILP repair.

[0074] Compared with directly solving the complete MILP, this method can reduce invalid integer searches and improve the solution efficiency and stability of large-scale petrochemical plant scheduling problems.

[0075] The abbreviations of the main production equipment in the text are shown in Table 1.

[0076] Table 1

[0077] This embodiment provides a petrochemical production scheduling optimization method based on decomposition coordination and topology analysis. Addressing the problems of existing centralized mixed-integer programming models—such as large numbers of variables, numerous integer decisions, and long direct solution times for the complete model under conditions of multiple units, multiple materials, and fine time granularity—this application divides the complete scheduling problem into multiple production unit sub-problems and blending and delivery sub-problems based on petrochemical production topology. This decomposition allows multiple scheduling sub-problems to be solved independently or in parallel, reducing the problem size of a single solution and providing a structural foundation for subsequent ADMM decomposition and coordination. To address the issue of inconsistent coupling relationships between sub-problems after production topology decomposition, such as intermediate material connections and component oil connections, this application employs ADMM to coordinate and iterate the production unit sub-problems and the blending and delivery sub-problems. In each iteration, each scheduling subproblem is solved separately. The consensus objective is updated and the coupling residual is calculated based on the values ​​of the coupling variables of each subproblem. The multiplier variables and penalty parameters are further updated to improve the consistency of cross-subproblem coupling relationships such as materials, inventory, and operating modes. To address the problem that decomposition coordination cannot directly guarantee a complete feasible solution in mixed integer scheduling problems, this application uses ADMM decomposition coordination as a heuristic coordination framework. Instead of directly piecing together the results of decomposed subproblems into the final solution, it identifies a stable integer skeleton from the results of multiple rounds of coordination iterations. Subsequently, based on the stable integer skeleton, constrained feasibility repair is performed in the complete mixed integer programming model, and the unfixed variables are re-optimized so that the final output solution meets the complete production constraints.

[0078] The method of the present invention will be further described below with reference to a specific embodiment. To verify the feasibility and effectiveness of the method described in the present invention, simulation experiments were conducted for three short-term oil refining production scheduling cases. Table 2 shows the basic information of the three test cases, which correspond to different order structures, scheduling loads, and logistics coupling complexities, respectively. The experiment used 1 hour as the main time step.

[0079] This time granularity can describe the process of plant operation, material transportation, inventory changes, product blending and order delivery in more detail. At the same time, it will increase the number of variables, binary variables and constraints in the model, thus reflecting the difficulty of solving large-scale petrochemical plant scheduling problems.

[0080] At a time step of 1 hour, the method of the present invention is compared with the global solution method of the solver, the Lagrange decomposition method and the rolling time domain decomposition method to verify the applicability of the method of the present invention under different scheduling cycles, order structures and logistics coupling complexities.

[0081] Table 2

[0082] To illustrate the model scale of the test cases, the number of variables and constraints for each test case under a 1-hour time step was further analyzed, and the results are shown in Table 3. The total number of variables for the three test cases is approximately 40,000 to 60,000, and the total number of constraints is approximately 90,000 to 120,000, with the number of binary variables all reaching tens of thousands.

[0083] Therefore, the test case belongs to the refining production scheduling problem with characteristics of large-scale mixed integer linear programming.

[0084] Table 3

[0085] The experiment uses the objective function value, solution time, and relative deviation as evaluation indicators, and the meanings of each indicator are shown in Table 4. The objective function value characterizes the overall cost of the scheduling scheme, including production and operating costs, inventory costs, and stockout penalties; the solution time characterizes the computation time required for the corresponding method to obtain a feasible scheduling scheme; and the relative deviation characterizes the difference between the objective function value obtained by the corresponding method and the optimal objective function value in this comparative experiment. For the scheduling problem described in this embodiment, the smaller the objective function value, the better the overall economic efficiency of the scheduling scheme; the shorter the solution time, the higher the solution efficiency of the corresponding method; and the smaller the relative deviation, the closer the objective function value obtained by the corresponding method is to the optimal objective function value in this comparative experiment.

[0086] Tables 4 to 6 show the solution results of the four methods in three test cases with a time step of 1 hour. The experimental results show that the method of this invention obtained the optimal objective function value in all three test cases in this comparative experiment, and the solution time was significantly reduced. This indicates that the method of this invention can generate feasible scheduling schemes with good overall economic efficiency under different scheduling scales and order structures.

[0087] The global solver method and the Lagrange decomposition method reached near the preset time limit of 3600 seconds in multiple test cases, indicating that at fine time granularity, the global solution of the complete MILP model and the traditional Lagrange decomposition solution still face significant computational pressure, and the obtained solution could not be proven to be the globally optimal solution within the time limit. The rolling time domain decomposition method has a shorter solution time, but its objective function value is relatively high, indicating that segmented optimization according to local time windows may reduce the coordination between consecutive time periods and affect the overall economy of the scheduling scheme.

[0088] In comparison, the method of this invention reduces the solution scale of individual subproblems by production topology decomposition, uses ADMM to coordinate the coupling relationship between materials, inventory and operation modes between devices, and uses restricted MILP repair to ensure that the final scheduling scheme meets the constraints of the complete model, thereby achieving a better balance between solution time and objective function value.

[0089] Table 4

[0090] Table 5

[0091] Table 6

[0092] The experimental results of the three test cases show that the method of the present invention is not only applicable to a single scheduling case, but also achieves the optimal objective function value in this comparative experiment under different scheduling periods, order quantities and delivery window structures.

[0093] Taking Case 3 as an example, this case has a long scheduling cycle, overlapping order delivery windows, and significant changes in demand for different products, reflecting the difficulty of production scheduling under complex coupling conditions. In this case, the solution time of the global solver method and the Lagrange decomposition method both reached or approached the preset time limit; although the rolling time domain decomposition method had a shorter solution time, the obtained objective function value was relatively high; the method of this invention obtained the feasible scheduling scheme with the best objective function value in the comparative experiment within 392.4 seconds. Experimental results show that under complex conditions of order delivery, inventory changes, and inter-unit logistics connections, the method of this invention can still obtain a scheduling scheme with good overall economic efficiency by utilizing the ADMM coordination process and the stable integer skeleton repair process.

[0094] At the scheduling scheme level, compared with the rolling time-domain decomposition method, the scheduling scheme obtained by the method of this invention has better continuity in terms of the switching of the operating mode of the main device and the adjustment of traffic flow, and can effectively avoid the mismatch problem between the preceding and following stages caused by independent optimization of local time windows. Compared with the global solver method and the Lagrange decomposition method, the method of this invention can obtain a feasible scheduling scheme that strictly satisfies the constraints of the complete model in a shorter computation time. The above results show that the combination of ADMM decomposition coordination strategy and stable integer skeleton repair mechanism is beneficial to improving the engineering executability of the scheduling scheme.

[0095] In addition to the main experiment with a time step of 1 hour, supplementary experiments with time steps of 2 hours, 4 hours, and 8 hours were also conducted to verify the applicability of the method of the present invention at different time granularities. Table 8 shows the solution results of the four methods at different time steps, using Case 3 as an example.

[0096] As shown in Table 7, the method of the present invention obtained the optimal objective function value in this comparative experiment under time steps of 2 hours and 4 hours, and the solution time was 72.4 seconds and 24.6 seconds, respectively. This indicates that the method of the present invention is not only applicable to fine-grained scheduling scenarios with a time step of 1 hour, but also to scheduling scenarios with a medium time step.

[0097] When the time step increases to 8 hours, the size of the complete scheduling model decreases accordingly, reducing the computational burden on the overall solver method and the traditional decomposition method. At this time step, the objective function value obtained by the method of this invention is not the optimal objective function value in this comparative experiment, and its advantage over other methods is weakened. These results indicate that the applicable advantages of the method of this invention are mainly reflected in petrochemical plant scheduling scenarios with fine to medium time granularity, large model size, and strong coupling relationships.

[0098] Table 7

[0099] Based on the simulation results above, the feasibility and effectiveness of the method of this invention have been verified. First, the three test cases cover different order quantities, scheduling cycles, demand scales, and delivery window structures, reflecting the short-term production scheduling needs of refining under varying degrees of complexity. Second, statistical results of the model scale show that, with a 1-hour time step, each test case contains tens of thousands of variables, tens of thousands of binary variables, and nearly one hundred thousand constraints, exhibiting typical characteristics of large-scale MILP scheduling problems. Third, the main experimental results show that the method of this invention can obtain feasible scheduling schemes that satisfy the complete model constraints in all three test cases, and the obtained objective function values ​​are the optimal objective function values ​​in this comparative experiment; its solution time is significantly lower than that of the global solver method and the Lagrange decomposition method, which reach the preset time limit. The representative scheduling scheme diagram also shows that the scheduling scheme obtained by the method of this invention has good continuity in terms of unit operation mode, flow adjustment, and logistics connection. Furthermore, supplementary experiments at different time steps show that the method of the present invention has good applicability in scheduling scenarios with fine to medium time granularity; when the time step increases and the model size decreases, the advantages of the method of the present invention compared with other methods are weakened.

[0100] Therefore, the method of the present invention can handle the inter-unit coupling relationship and a large number of integer decisions in the scheduling of large-scale petrochemical units. On the basis of meeting production constraints such as material balance, inventory changes, unit operation, product blending and order delivery, it can obtain a scheduling scheme with good comprehensive economics and engineering feasibility, and has the feasibility to be applied to petrochemical production scheduling scenarios.

[0101] It should be understood that various parts of the present invention can be implemented in hardware, software, firmware, or a combination thereof.

[0102] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A petrochemical production scheduling optimization method based on decomposition coordination and topology analysis, characterized in that, include: Step 1: Obtain the production units, discrete time periods, operating modes, material flow, inventory, product blending, and order delivery data of the petrochemical production system to be scheduled, form a production scheduling instance, and establish a production topology based on the connection relationships between production units, materials, virtual storage tanks, product blending links, and order delivery links. Step 2: Based on the production scheduling example, construct a discrete-time global model for petrochemical production scheduling. The global model aims to minimize the comprehensive cost within the scheduling cycle. At the production scheduling level, it describes the constraints between the operating status of production units, mode switching process, material input and output, inventory changes, product blending, and order delivery. The global model is a mixed-integer linear programming model. Step 3: Based on the production topology, the scheduling problem corresponding to the global model is divided into multiple production unit sub-problems and one blending and delivery sub-problem. Each production unit sub-problem corresponds to one production unit, and the blending and delivery sub-problem retains the component oil inventory, finished oil inventory, blending flow rate, blending ratio, quality specifications, and order delivery constraints. Step 4: Transform the cross-subsystem connection relationships in the global model into coordination residuals, wherein the coordination residuals include material-oil balance residuals, component-oil connection residuals, and coupling device mode residuals, which are used to quantify the degree of deviation of the solutions of each subproblem from the cross-boundary coupling relationships. Step 5: Perform linear relaxation initialization, relax the integer variables in the global model into continuous variables and solve the linear relaxation model to obtain a reference solution that reflects the global scheduling trend. After rounding the corresponding binary variables in the reference solution, perform the first round of feasibility repair on the integer trajectory of the fixed entry device to obtain the initial feasible upper bound. Step 6: Perform decomposition coordination iteration guided by the alternating direction multiplier method. In each round of iteration, solve the subproblems of each production unit and the coordination and delivery subproblems respectively. Calculate the original residual of this round according to the definition of coordination residual, update the coordination objective, penalty parameters and Lagrange multipliers, and perform multi-round iterative coordination. In step 6, the penalty parameter is adaptively updated based on the relative relationship between the original residual and the dual residual. Let the first The original residual norm and the dual residual norm are respectively and The residual ratio is defined as: ,when Located in the preset balance range Internal time maintenance ,in The imbalance threshold parameter; when When the value exceeds the preset balance range, the penalty parameter is updated as follows: in, To update the index proportionally, To smooth out the weights, and These are the limiting coefficients for single-wheel increase and single-wheel decrease, respectively. and These are the lower and upper bounds of the penalty parameter, respectively. Indicates will Restricted to closed intervals Inside; Step 7: During the iterative process of the alternating direction multiplier method, continuously record the integer pattern trajectories returned by the subproblems of each production unit, and generate candidate fixed skeletons based on the stability of the steady-state pattern variables of each unit in the most recent multiple iterations. Step 7 specifically includes: Suppose the system is in the most recent The trajectory of the steady-state variable recorded in the round of iteration is For the device , its first wheel stability Defined as: in, Indication device The set of steady-state variables, Indicates the index of the steady-state variable. The cardinality of a set is the number of elements contained in the set. This is an indicator function that takes the value 1 if the condition within the parentheses is true, and 0 otherwise. For devices with the same stability, the device trajectory is selected first to enter the candidate fixed frame in order of the original residual values ​​from low to high. Step 8: Construct a restricted mixed integer linear programming repair model based on the candidate fixed skeleton, fix the integer trajectories of some stable devices, and re-optimize the integer variables of the unfixed devices as well as all continuous logistics, inventory, re-mixing and delivery variables to obtain a feasible scheduling scheme that satisfies all constraints of the global model. Step 9: Repeat steps 6 to 8 until the termination condition is met, and output the final scheduling scheme.

2. The method according to claim 1, characterized in that, The production device subproblem retains the device's own steady-state mode variables, transition start-up variables, transition occupancy variables, feed boundary, yield relationship, and local production cost. The transition start-up variable is used to indicate whether the production unit has started the switching process from the source operating mode to the target operating mode, and the transition occupancy variable is used to indicate whether the production unit is in a transition occupancy state during the switching from the source operating mode to the target operating mode.

3. The method according to claim 1, characterized in that, The material oil balance residual The expression is: in, Indicates production material oil A collection of devices, Indicates the consumption of material oil A collection of devices, Indication device During the period Input flow, Indication device During the period For material oil ; output flow rate; The component oil connection residual The expression is: in, Indicates the production of component oils A collection of devices, This indicates the time period for reconciliation and delivery. For component oil Input volume, Indication device During the period Producing component oil Traffic, and These represent the component oils in the virtual tank. During the period and time period Inventory levels Indicates the blending and delivery side component oils The change in inventory of virtual containers between adjacent time periods; The coupling device mode residual The expression is: in, and Representing the devices and device During the period Is it in mode? .

4. The method according to claim 1, characterized in that, The coordination target is updated using a historical smoothing method, and the update formula is: in, Indicates the reconciliation target index. This indicates the updated coordination objective. Indicates the first The coordination goal of the wheel Indicates historical weight, Indicates by the first The original coordination value obtained from the solution to the wheel problem.

5. The method according to claim 1, characterized in that, The update method for the Lagrange multipliers is as follows: in, Indicates the first The coordination residuals that are due for them Indicates the first The Lagrange multipliers of the wheel, Indicates the first The Lagrange multipliers of the wheel have not yet undergone scalar adjustment; If the penalty parameter changes, the multiplier is adjusted for synchronization scaling: in, and They represent the first Round iteration and the first Wheel penalty parameters, This indicates the number obtained after scale adjustment. The Lagrange multipliers of the wheel; If the penalty parameter does not change, then .

6. The method according to claim 1, characterized in that, The processing steps for the restricted mixed-integer linear programming repair model include: If the current candidate skeleton is different from the previous skeleton, then perform a restricted mixed integer linear programming repair based on the current candidate skeleton, and compare the repaired objective function value with the historical best upper bound. If it is better than the historical best upper bound, then update the historical best feasible solution and the corresponding objective value. If the current candidate skeleton is exactly the same as the previous skeleton, skip this round of repair.

7. The method according to claim 1, characterized in that, The termination condition is: The historical optimal upper bound has not improved within several consecutive rounds, the integer trajectories of each device tend to stabilize, the candidate skeleton appears repeatedly, or the improvement of the original residual enters a plateau period, or the preset maximum number of iterations is reached.

Images