CS 30 Notes on Using Matrix Operations to Model the Flow of Materials in a Refinery

Mart Molle, April 2011

Background

This material is based on my experience as an undergraduate student intern during the summer of 1975, doing computer programming to support the engineers managing a copper refinery in northern Canada. I have forgotten plenty of details, and simplified many others to make the example easier to understand. Please don't try to construct a smelter based on these instructions! For the purposes of this example, let's assume that refining copper ore into finished copper bars requires two processing steps, as shown in the following flow diagram.
Flow Diagram for Copper Smelter

Step 1: Blast Furnace (simple diagram)

The smelting process begins by loading some raw materials (copper ore mixed with other inputs, such as silica and coal) into a blast furnace. Then lots of heast is applied until the mixture melts and then separates into several layers -- just like oil and water in a bottle of salad dressing, or the fat and meat juices you collect from the bottom of the roasting pan when making gravy. In this case, the top layer is impurities suspended in molten silica, called slag, which cools into a solid material that looks like lava and is discarded as waste. The bottom layer is rough copper, which is cast into slabs in preparation for the second step. In between is an undifferentiated mixture of ingredients with similar composition to the raw materials.

For the purpose of this example, let's assume that one "batch" of raw materials at a time is processed by the blast furnace:
In addition, let's assume that the refining process carried out by the furnace distributes the input copper among the three outputs independently of all other material present inside the smelter. In other words, looking only at the processing of copper in our example, we see:

Step 2: Electrolysis Tank

The second step of the smelting process involves electro-deposition, which is similar charging a battery or applying chrome plating to the exhaust pipes of a motorcycle. In this case, the rough copper bars are immersed in a tank of electrolyte next to some sheets of refined copper. When a DC voltage is applied between the anode (rough copper) and cathode (refined copper), the copper atoms move through the electrolyte as charged ions from anode to cathode, until the anode is reduced to some left-over sludge (containing 10% copper, 90% other) sitting on the bottom of the tank and the cathode has grown into a large piece of refined copper (containing 99% copper, 1% other).

Note, that this information only tells us the relative amount of copper in each output product, without giving you the total weight of each output the way we had for the blast furnace. However, we can find this missing incredient by solving a "mass balance" equation for tracking the movement of copper. Assuming that no material is created or destroyed inside the electolysis tank, our input of 10 tons of raw copper must generate a total of 10 tons of output, which is split into z and (10 - z) tons, respectively, of sludge and refined copper. Similarly, the total amount of copper (in tons) entering and leaving the electrolysis tank is also equal, which gives us the following equation for the movement of copper:

10 * 0.75 = z * 0.1 + (10 - z) * 0.99

This equation can be easily solved algebraically to give

z = 2.4 / 0.89 = 2.697

(Remember that the above lines are equations, not Matlab commands, and the first line in particular will generate a Matlab error, rather than doing what you wanted!)

Therefore, once we have found z, we can represent the processing of copper at the electrolysis tank as follows:
In other words, the processing in the electrolysis step distributes the input copper so that

(0.27 tons) / (7.5 tons) = 3.6%

is caught in the sludge, and the remaining

(7.23 tons) / (7.5 tons) = 96.4%

leaves the smelter in the form of its primary product: refined copper bars.

It is also very important to not to forget about the sludge, which is an "unwanted" output product generated by the electrolysis tank. Since the sludge the same concentration of copper (10%) as the original raw materials entering the smelter, the sludge must be returned to the blast furnace for further processing, rather than simply thrown away.

Modeling the Smelter with Matlab

After a few preliminary calculations (i.e., solving for z), we have converted the operation of the smelter into a series of simple formulas, to say how the amount of copper (say) in each input/output stream or processing step depends on the others. Thus, we could easily write down a series of Matlab assignment statements to show how these values change from one shift to the next. For example, suppose the smelter is initially idle after a long vacation, and the first shift shows up for work and starts to bring the smelter back online. Since both the blast furnance and electrolysis tank are empty, the only crew members at the input stockpile can do any real work, namely moving 10 tons of raw materials to the entrance of the blast furnace. Thus, a "snapshot" of the smelter taken during the first shift can be summarized by the following Matlab statement:

% first shift
input_copper = 10;    % prepare first batch of raw material

Similarly a smelter "snapshot" during the second shift looks like this. The blast furnace crew has loaded the input copper into the furnace and is now heating it.  The stockpile crew is moving another 10 tons of raw materials to the entrance of the blast furnace. All other crew members are still idle, so the new "state" of the smelter looks like this:

% second shift
blast_furnace = input_copper
;    % load blast furnace
input_copper = 10;                        % prepare another batch of raw material
 
As the process continues, we get the follow sequence of Matlab statements:

% third shift
slag = blast_furnace * 0.0;    % empty the blast furnace
raw_copper = blast_furnace * 0.75;
leftovers = blast furnace * 0.25;

blast_furnace = leftovers + input_copper; % reload the blast furnace
input_copper = 10;    % prepare another back of raw material

tank = raw_copper;    % load the electrolysis tank

slag_dump = slag;    % disposal crew takes slag to dump
slag = 0;

% fourth shift
slag = blast_furnace * 0.0;    % empty the blast furnace
raw_copper = blast_furnace * 0.75;
leftovers = blast furnace * 0.25;

sludge = raw_copper * 0.036;    % empty the electrolysis tank
refined_copper = tank * 0.964;

tank = raw_copper;    % reload the electrolysis tank
blast_furnace = leftovers + input_copper + sludge; % reload the blast furnace
input_copper = 10;    % prepare another back of raw material

output_product = refined_copper;    % shipping crew sends final product to market
refined_copper = 0

slag_dump = slag;    % disposal crew takes slag to dump
slag = 0;

% fifth shift ...  Looks exactly the same as the fourth shift

Simplifying the Model by Introducing Matlab Matrix/Vector Notation.

Vector notation allows us to combine the results of one processing "step" at the blast furnace, i.e.,

slag = blast_furnace * 0.0;    % empty the blast furnace
raw_copper = blast_furnace * 0.75;
leftovers = blast furnace * 0.25;

into a single Matlab command:

blast_out = blast_furnace * [0.0 0.75 0.25]

blast_out =    0   7.5000   2.5000

(where I have shown the answer produced by Matlab in red to separate it from the commands I've been typing). The only "inconvenience" is that I must recognize each output product by its position in the vector, rather than a "nice" individually-chosen variable. Don't worry, things will make more sense in a moment!

Notice that all the objects in the smelter diagram has been numbered from 1 to 5. Therefore, let's begin by using those numbers to identify each smelter object by its position in a 5-element, as shown in the following table:

Old Name:
input_copper
blast_furnace
slag_dump
tank
output_product
New Name:
smelter(1)
smelter(2)
smelter(3)
smelter(4)
smelter(5)

For example, the current state of the smelter might look like this after running for a few shifts

smelter

smelter =
   10.0000   13.3950         0    9.3750    7.2300

Furthermore, notice that each of the "streams" of material (raw input, slag, etc) has a comes from a single "producer" and goes to a single "consumer". Thus, we can identify each "stream" produced by an object by the number assigned to its consumer. In other words, even though the blast furnace produces only three outputs, we should call them blast_out(1) through blast_out(5) so we know where they are going.  In this case, our previous Matlab command would be rewritten like this:

blast_out = smelter(2) * [0.0 0.25 0.0 0.75 0.0]

blast_out =
        0   2.5000   0
   7.5000  0

Similarly, we would represent the result of one processing step at the electrolysis tank as:

tank_out = smelter(4) * [0.0 0.036 0.0 0.0 0.964]

tank_out =
        0   0.2700   0
  0   7.2300

For the remaining locations in the smelter, the result of one processing step is very easy. Everything at the input stockpile is loaded into the blast furnce, like this:

stockpile_out = smelter(1) * [0.0 1.0 0.0 0.0 0.0]

stockpile_out =
        0   10.0000   0   0   0


(I know the workers will refill the stockpile with a new input material, but right now I'm focussing my attention on the movment of existing material within the smelter. We will return come back to the input in a few minutes.) On the other hand, everything at the slag dump and product shipping area is permanently removed from the system. Since none of the material currently at either smelter(3) and smelter(5) is still somewhere inside the smelter one shift from now, we can represent the those parts of the process as a multiplication by a vector of zeros.

Let us now stack up these "processing vectors" for each smelter location to form a two dimensional matrix with location 1 (input stockpile) forming row 1 and the product shipping area forming row 5:

shift =  [0    1.0    0    0    0
              0    0.25    0    0.75    0
              0    0    0    0    0
              0    0.036    0    0    0.964
              0    0    0    0    0]

Notice that we can use the colon operator to select a specific row from this matrix:

shift(4,:)

ans =
         0    0.0360         0         0    0.9640

and thus we can rewrite the previous command for the electrolysis tank as:

tank_out = smelter(4) * shift(4,:)

tank_out =    0   0.2700   0  0   7.2300

More generally, let us now repeat this "trick" to generate the outputs from one processing step at each of the five objects in the smelter model:

stockpile_out = smelter(1) * shift(1,:);    % notice the semicolons to suppress output
blast = smelter(2) * shift(2,:);
slag_out = smelter(3) * shift(3,:);
tank_out = smelter(4) * shift(4,:);
product_out = smelter(5) * shift(5,:);

Clearly, we can find the input to some "target" object at the next processing step, by summing all outputs generated in the current processing step that are sent directly to the "target" object. For example, if object 2 (the blast furnace) is our "target" object, then its input for the next shift must be:

stockpile_out(2) + blast(2) + slag_out(2) + tank_out(2) + product_out(2)

Fortunately, we can use matrix operations to avoid having to deal with all these tedious details. More specifically, let's now go back and look at the outputs generated by the "trick" a few lines ago. (Notice how I used the disp( ) command, so we see only the values of each variable, which is a  5-element row vector, without showing any of the variable names.)

disp(stockpile_out), disp(blast), disp(slag_out), disp(tank_out), disp(product_out)

     0    10    0     0     0

     0    3.3487    0   10.0463   0

     0     0     0     0     0

     0    0.3375    0     0    9.0375

     0     0     0     0     0

To convert this two-dimensional table of values into the contents of the kth object in the smelter model at the next processing step, we must add together every element in the kth column. Using vector notation, we can do this simultaneously for every object in a single Matlab expression, like this:

stockpile_out + blast + slag_out + tank_out + product_out
ans =
   10.0000   13.6862         0   10.0463    9.0375



So here is the bottom line for all of the discussion above. The way we got from the 5-element row vector smelter and  5x5 matrix shift was to multiply the rows of shift by the corresponding element from the vector smelter to form a 5x5 table of results, and then "flattening" the table by adding the elements in each column. This entire series of steps is exactly what happens automatically when we apply ordinary matrix multiplication to vector smelter and array shift.

Thus, all we really needed to do to model the operation of the smelter is the following. To move forward in time by one shift, we must redistribute the contents of smelter (i.e., the material currently in the system) by multiplication by the shift matrix, then adding the the vector input (representing the new material added from the stockpile every shift), like this:

>> input=[10 0 0 0 0], smelter=input

input =

    10     0     0     0     0


smelter =

    10     0     0     0     0

>> smelter=smelter*shift+input

smelter =

    10    10     0     0     0

>> smelter=smelter*shift+input

smelter =

   10.0000   12.5000         0    7.5000         0

>> smelter=smelter*shift+input

smelter =

   10.0000   13.3950         0    9.3750    7.2300

>> smelter=smelter*shift+input

smelter =

   10.0000   13.6862         0   10.0463    9.0375

>> smelter=smelter*shift+input

smelter =

   10.0000   13.7832         0   10.2647    9.6846

>>

We can keep going forward like this, one shift at a time, way for as long as we like. However, the obvious question is what happens a very long time in the future? Is there a long-term "steady-state" distribution of material in the smelter over a shift, and if so, how do we find it?

The answer to this question is very simple in Matlab, if you apply what you know about how to find the sum of a (scalar) geometric series.  In other words, suppose I make n new friends every day, but at the same time loose a fraction (1-p) of my existing friends because I don't pay enough attention to them. How many total friends will I have in the distant future?  Clearly I will have n from today, n*p from yesterday, n*p^2 from two days ago, etc, giving me a total of n/(1-p). Similarly, the contents of the smelter will include input from the current shift, input*shift remaining from the previous shift, (input*shift)*shift remaining from two shifts ago, etc., which is just another geometric series, which can be solved in the same way to give:

>> input/(eye(5)-shift)

ans =

   10.0000   13.8313         0   10.3734   10.0000

>>