• Source Code for Ternary Object Module and Combinatorial Analysis Procedures
• Enter Ternary System Data
• Plot the Ternary System
• Flowsheet Synthesis
  • Enter composition points and units
  • Solve for Feasible Flowsheets
  • View Flowsheet(s)
• Restart

Synthesis of Separation Systems for Ternary Mixtures

Lee R. Partin

Copyright 2000-2001, L R Partin Enterprises

September 8, 2001

lpartin@chartertn.net

http://users.chartertn.net/lpartin

This application combines the power of residue curve maps (RCMs) and combinatorial analysis for the synthesis of separation systems for ternary mixtures.  RCMs provide knowledge of feasible distillation separations.  The skilled engineer can interpret the RCM to define distillation towers and related unit operations for separating azeotropic mixtures.  Combinatorial analysis is combined with the unit operation definitions to handle the synthesis task.  The MSG algorithm determines when a feasible flowsheet superstructure (maximal structure) has been established.  The SSG algorithm then extracts all of the feasible flowsheets from the maximal structure.

Source Code for Ternary Object Module and Combinatorial Analysis Procedures

>

with(plottools): with(networks): with(combinat): # add a special triangular map coordinate system addcoords(triangular,[x,y],[x+y*cos(Pi/3),y*cos(Pi/6)]): # create a ternary map object module T()   description "Ternary DBM and Flowsheeting Object";      local i, Sequences, TLval, TIval, THval, TLIval, TIHval, TLHval, TLIHval, MapSelection, TSet, p1, p2, p3, p4, p5, p6, baseplot, baseplot2, AzLIval, AzIHval, AzLHval, AzLIHvalL, AzLIHvalI, MoveToTriangular, plotNodeSaddle, plotBoundaries, Boundaries, points, plotPoints, Psi_minus, Psi_plus, phi_minus, phi_plus,  SSGmethod,  Raw_Materials, Products, by_prods, Op, decision_map,  ox_vals, Units, maximal_structure, Pgraph;      export  TSequences, DBM, TL, TI, TH, TLI, TIH, TLH, TLIH, checkT, printData,  TSequence, PlotDBM, AzLI, AzIH, AzLH, AzLIH, MapSelectionNumber, AddPoint, PrintPoints, AddUnit, MSG, SSG, graph_network, flowsheet, PlotPoints, NumberOfMaps, MapUnits;     # the map selection is initially set to 1.  when there are more than   # one feasible map for a temperature sequence, then the first map   # is used.   MapSelection:=1;   # initialize the position temperatures   TLval:=NULL;   TIval:=NULL;   THval:=NULL;   TLIval:=NULL;   TIHval:=NULL;   TLHval:=NULL;   TLIHval:=NULL;   NumberOfMaps:=NULL;   MapUnits:="Units: ";   TSet:=[TLval,TIval,THval];   points:={}; Raw_Materials:={}; Products:={}; Units:={};   decision_map:={}; maximal_structure:={};      # Boundaries is a table of DBM boundary lines Boundaries:=table(["000"=[],"001"=[],"002-m"=[67],"003"=[],"004-M"=[67],"010-S"=[17,37,47,57],"011-S"=[37,47,57,67],"012"=[46],"013-S"=[17,37,47,67],"020"=[14],"021"=[46],"022-m"=[47,67],"023"=[14],"024-M"=[14,67],"030"=[],"031"=[],"032-m"=[67],"033-S"=[17,47,57,67],"034"=[46],"040-M"=[47],"041-M"=[47],"043"=[46],"044-M"=[47,67],"100"=[],"101-S"=[17,27,57,67],"102"=[26],"103"=[],"104-M"=[67],"110-S"=[27,37,47,57],"112-S"=[27,37,47,67],"112-SH"=[26,27,37,47,57],"112-SL"=[17,27,37,47,46],"113-S"=[27,37,47,67],"120"=[24],"121-S"=[17,27,47,67],"121-SH"=[17,27,24,57,67],"121-SI"=[17,27,37,46,67],"122"=[24,26],"123"=[24],"124-M"=[24,67],"130"=[],"131-S"=[17,27,47,67],"132"=[26],"133-S"=[27,47,57,67],"134"=[46],"140-M"=[47],"142-M"=[26,47],"143"=[46],"144-M"=[47,67],"200-m"=[27],"201"=[26],"202-m"=[27,67],"203-m"=[27],"210"=[24],"211-S"=[27,47,57,67],"211-SI"=[26,37,47,57,67],"211-SL"=[17,24,47,57,67],"212"=[24,46],"213"=[24],"214-M"=[24,67],"220-m"=[27,47],"221"=[26,46],"222-m"=[27,47,67],"223-m"=[27,47],"230-m"=[27],"231"=[26],"232-m"=[27,67],"234-m"=[27,46],"241-M"=[26,47],"243-m"=[27,46],"300-S"=[17,27,37,57],"301-S"=[27,37,57,67],"303-S"=[17,27,37,67],"304"=[26],"310-S"=[17,27,47,57],"311-S"=[27,47,57,67],"312"=[46],"313-S"=[17,27,47,67],"314"=[26],"320"=[14],"321"=[46],"322-m"=[47,67],"323"=[14],"324-M"=[14,67],"324-m"=[26,47],"330-S"=[17,27,37,47],"331-S"=[27,37,47,67],"334-S"=[27,37,47,67],"334-SH"=[26,27,37,47,57],"334-SL"=[17,27,37,47,46],"340"=[24],"341"=[24],"342-m"=[24,67],"343-S"=[17,27,47,67],"343-SH"=[17,27,24,57,67],"343-SI"=[17,27,37,46,67],"344"=[24,26],"400"=[25],"401"=[25],"402-m"=[25,67],"403"=[26],"404-M"=[27,67],"410"=[25],"411"=[25],"412-M"=[27,46],"412-m"=[25,67],"413"=[26],"414-M"=[27,67],"420-M"=[14,27],"420-m"=[25,47],"421-M"=[27,46],"421-m"=[25,47],"423-M"=[14,27],"423-m"=[26,47],"430"=[24],"431"=[24],"432-m"=[24,67],"433-S"=[27,47,57,67],"433-SI"=[26,37,47,57,67],"433-SL"=[17,24,47,57,67],"434"=[24,46],"440-M"=[27,47],"441-M"=[27,47],"443"=[26,46],"444-M"=[27,47,67]]);   # Sequences stores the DBM names for each valid temperature sequence   Sequences:=[[135, `000`], [1325, `400`], [1354, `030`], [1356, `003`],   [1435, `020`], [2135, `100`], [6135, `001`], [13254, `430`], [13256,   `403`], [13524, `430`], [13526, `403`], [13542, `340`], [13546, `043`],   [13547, `040-M`], [13562, `304`], [13564, `034`], [13567, `004-M`],   [13725, `300-S`], [13752, `300-S`], [14325, `320 410 `], [14352, `320`], [14356, `023`], [14735, `010-S`], [21354, `130`], [21356, `103`], [21435, `120`], [24135, `120`], [26135, `102`], [41325, `410`], [41735, `010-S`], [42135, `210`], [46135, `012`], [61325, `401`], [61354, `031`], [61435, `021`], [62135, `201`], [64135, `021`], [72135, `200-m`], [76135, `002-m`], [132546, `443`], [132547, `440-M`], [132564, `434`], [132567, `404-M`], [135246, `443`], [135247, `440-M`], [135264, `434`], [135267, `404-M`], [135426, `443`], [135427, `440-M`], [135462, `344`], [135467, `044-M`], [135624, `434`], [135627, `404-M`], [135642, `344`], [135647, `044-M`], [135724, `330-S`], [135726, `303-S`], [135742, `330-S`], [135746, `033-S`], [135762, `303-S`], [135764, `033-S`], [137254, `330-S`], [137256, `303-S`], [137524, `330-S`], [137526, `303-S`], [137542, `330-S`], [137562, `303-S`], [143256, `323 413 `], [143257, `420-M`], [143275, `420-M`], [143526, `323 413 `], [143527, `420-M`], [143562, `314 323 `], [143567, `024-M`], [143725, `310-S`], [143752, `310-S`], [147325, `310-S`], [147352, `310-S`], [147356, `013-S`], [174325, `420-m`], [213546, `143`], [213547, `140-M`], [213564, `134`], [213567, `104-M`], [214356, `123`], [214735, `110-S`], [241356, `123`], [241735, `110-S`], [246135, `122`], [247135, `110-S`], [261354, `132`], [261435, `122`], [264135, `122`], [267135, `101-S`], [413256, `413`], [413526, `413`], [413562, `314`], [413725, `310-S`], [413752, `310-S`], [417325, `310-S`], [417352, `310-S`], [417356, `013-S`], [421356, `213`], [421735, `110-S`], [426135, `212`], [427135, `110-S`], [461325, `312 411 `], [461352, `312`], [461735, `011-S`], [462135, `212`], [467135, `011-S`], [613254, `431`], [613524, `431`], [613542, `341`], [613547, `041-M`], [613725, `301-S`], [613752, `301-S`], [614325, `321 411 `], [614352, `321`], [614735, `011-S`], [621354, `231`], [621435, `221`], [624135, `221`], [627135, `101-S`], [641325, `321 411 `], [641352, `321`], [641735, `011-S`], [642135, `221`], [647135, `011-S`], [714325, `420-m`], [721354, `230-m`], [721356, `203-m`], [721435, `220-m`], [724135, `220-m`], [726135, `202-m`], [741325, `420-m`], [742135, `220-m`], [746135, `022-m`], [761325, `402-m`], [761354, `032-m`], [761435, `022-m`], [762135, `202-m`], [764135, `022-m`], [1325467, `444-M`], [1325647, `444-M`], [1325746, `433-S 433-SI 433-SL `], [1325764, `433-S 433-SI 433-SL `], [1352467, `444-M`], [1352647, `444-M`], [1352746, `433-S 433-SI 433-SL `], [1352764, `433-S 433-SI 433-SL `], [1354267, `444-M`], [1354627, `444-M`], [1354726, `343-S 343-SH 343-SI `], [1354762, `343-S 343-SH 343-SI `], [1356247, `444-M`], [1356427, `444-M`], [1356724, `334-S 334-SH 334-SL `], [1356742, `334-S 334-SH 334-SL `], [1357246, `343-SI 433-SI 433-SL `], [1357264, `334-SL 433-SI 433-SL `], [1357426, `343-SH 343-SI 433-SI `], [1357462, `334-SH 343-SH 343-SI `], [1357624, `334-SH 334-SL 433-SL `], [1357642, `334-SH 334-SL 343-SH `], [1372546, `343-SI `], [1372564, `334-SL `], [1375246, `343-SI `], [1375264, `334-SL `], [1375426, `343-SI `], [1375462, `343-SI `], [1375624, `334-SL `], [1375642, `334-SL `], [1432567, `324-M 414-M 423-M `], [1432576, `423-M`], [1432756, `423-M`], [1435267, `324-M 414-M 423-M `], [1435276, `423-M`], [1435627, `324-M 414-M 423-M `], [1435672, `324-M`], [1435726, `313-S`], [1435762, `313-S`], [1437256, `313-S`], [1437526, `313-S`], [1437562, `313-S`], [1473256, `313-S`], [1473526, `313-S`], [1473562, `313-S`], [1743256, `423-m`], [1743526, `423-m`], [1743562, `324-m`], [2135467, `144-M`], [2135647, `144-M`], [2135746, `133-S`], [2135764, `133-S`], [2143567, `124-M`], [2147356, `113-S`], [2413567, `124-M`], [2417356, `113-S`], [2461735, `112-SH `], [2467135, `112-SH 112-SL 121-SH `], [2471356, `113-S`], [2476135, `112-S 112-SH 112-SL `], [2613547, `142-M`], [2614735, `112-SH `], [2641735, `112-SH `], [2647135, `112-SH 121-SH 121-SI `], [2671354, `131-S`], [2671435, `121-S 121-SH 121-SI `], [2674135, `121-S 121-SH 121-SI `], [4132567, `414-M`], [4135267, `414-M`], [4135627, `414-M`], [4135726, `313-S`], [4135762, `313-S`], [4137256, `313-S`], [4137526, `313-S`], [4137562, `313-S`], [4173256, `313-S`], [4173526, `313-S`], [4173562, `313-S`], [4213567, `214-M`], [4217356, `113-S`], [4261735, `112-SH `], [4267135, `112-SH 112-SL 211-SL `], [4271356, `113-S`], [4276135, `112-S 112-SH 112-SL `], [4613257, `412-M`], [4613275, `412-M`], [4613527, `412-M`], [4613725, `311-S`], [4613752, `311-S`], [4617325, `311-S`], [4617352, `311-S`], [4621735, `211-SI `], [4627135, `112-SL 211-SI 211-SL `], [4671325, `311-S`], [4671352, `311-S`], [4672135, `211-S 211-SI 211-SL `], [4761325, `412-m`], [6132547, `441-M`], [6135247, `441-M`], [6135427, `441-M`], [6135724, `331-S`], [6135742, `331-S`], [6137254, `331-S`], [6137524, `331-S`], [6137542, `331-S`], [6143257, `421-M`], [6143275, `421-M`], [6143527, `421-M`], [6143725, `311-S`], [6143752, `311-S`], [6147325, `311-S`], [6147352, `311-S`], [6174325, `421-m`], [6213547, `241-M`], [6214735, `211-SI `], [6241735, `211-SI `], [6247135, `121-SH 121-SI 211-SI `], [6271354, `131-S`], [6271435, `121-S 121-SH 121-SI `], [6274135, `121-S 121-SH 121-SI `], [6413257, `421-M`], [6413275, `421-M`], [6413527, `421-M`], [6413725, `311-S`], [6413752, `311-S`], [6417325, `311-S`], [6417352, `311-S`], [6421735, `211-SI `], [6427135, `121-SI 211-SI 211-SL `], [6471325, `311-S`], [6471352, `311-S`], [6472135, `211-S 211-SI 211-SL `], [6714325, `421-m`], [6741325, `421-m`], [7143256, `423-m`], [7143526, `423-m`], [7143562, `324-m`], [7213546, `243-m`], [7213564, `234-m`], [7214356, `223-m`], [7241356, `223-m`], [7246135, `222-m`], [7261354, `232-m`], [7261435, `222-m`], [7264135, `222-m`], [7413256, `423-m`], [7413526, `423-m`], [7413562, `324-m`], [7421356, `223-m`], [7426135, `222-m`], [7461325, `322-m 412-m 421-m `], [7461352, `322-m`], [7462135, `222-m`], [7613254, `432-m`], [7613524, `432-m`], [7613542, `342-m`], [7614325, `322-m 412-m 421-m `], [7614352, `322-m`], [7621354, `232-m`], [7621435, `222-m`], [7624135, `222-m`], [7641325, `322-m 412-m 421-m `], [7641352, `322-m`], [7642135, `222-m`]];          # p1 and p2 are prepared plots for ternary diagrams.   baseplot:={[[0,0],[1,0]],[[0,0],[0,1]],[[1,0],[0,1]]}: p1:=plot(baseplot,axes=`NONE`,color=`black`,coords=`triangular`,scaling=`CONSTRAINED`):  baseplot2:={seq([[i/10,0],[i/10,1-i/10]],i=1..9),seq([[0,i/10],[1-i/10,i/10]],i=1..9),seq([[i/10,0],[0,i/10]],i=1..9)}:   p2:=plot(baseplot2,axes=`NONE`,color=`black`,coords=`triangular`,linestyle=3,color=grey,scaling=`CONSTRAINED`):      # the list is converted into a table on sequence numbers.   TSequences:=table();   for i from 1 to nops(Sequences) do TSequences[Sequences[i][1]]:=       op(sscanf(Sequences[i][2],"%s %s %s")) end do;   TL:=proc() description "update low boiler T";              if nargs > 0 then                 if type(args[1],`numeric`) then TLval:=args[1]                 else error "Expected a numeric value,%1",args[1] end if              else TLval:=NULL end if;              TSet:=[TLval,TIval,THval,TLIval,TIHval,TLHval,TLIHval];              if TLval <> NULL then printf("TL = %g",TLval) end if;              checkT();   end proc;   TI:=proc() description "update intermediate boiler T";              if nargs > 0 then                 if type(args[1],`numeric`) then TIval:=args[1]                 else error "Expected a numeric value,%1",args[1] end if              else TIval:=NULL end if;              TSet:=[TLval,TIval,THval,TLIval,TIHval,TLHval,TLIHval];              if TIval <> NULL then printf("TI = %g",TIval) end if;              checkT();   end proc;     TH:=proc() description "update heavy boiler T";              if nargs > 0 then                 if type(args[1],`numeric`) then THval:=args[1]                 else error "Expected a numeric value,%1",args[1] end if              else THval:=NULL end if;              TSet:=[TLval,TIval,THval,TLIval,TIHval,TLHval,TLIHval];              if THval <> NULL then printf("TH = %g",THval) end if;              checkT();   end proc;   TLI:=proc() description "update LI azeotrope";              if nargs > 0 then                 if type(args[1],`numeric`) then TLIval:=args[1]                 else error "Expected a numeric value,%1",args[1] end if              else TLIval:=NULL end if;              if TLval<>NULL and TIval<>NULL and TLIval<>NULL and TLIval>TLval and TLIval<TIval then error "Bad TLI entry:  TLI must not be between TL and TI." end if;              TSet:=[TLval,TIval,THval,TLIval,TIHval,TLHval,TLIHval];              if TLIval <> NULL then printf("TLI = %g",TLIval) end if;              if nops(TSet)>nops({op(TSet)}) then                 error "The temperatures must be unique.  Two values are the same."                end if;   end proc;   TIH:=proc() description "update IH azeotrope";              if nargs > 0 then                 if type(args[1],`numeric`) then TIHval:=args[1]                 else error "Expected a numeric value,%1",args[1] end if              else TIHval:=NULL end if;              TSet:=[TLval,TIval,THval,TLIval,TIHval,TLHval,TLIHval];              if TIHval <> NULL then printf("TIH = %g",TIHval) end if;              if nops(TSet)>nops({op(TSet)}) then                 error "The temperatures must be unique.  Two values are the same."                end if;   end proc;     TLH:=proc() description "update LH azeotrope";              if nargs > 0 then                 if type(args[1],`numeric`) then TLHval:=args[1]                 else error "Expected a numeric value,%1",args[1] end if              else TLHval:=NULL end if;              TSet:=[TLval,TIval,THval,TLIval,TIHval,TLHval,TLIHval];              if TLHval <> NULL then printf("TLH = %g",TLHval) end if;              if nops(TSet)>nops({op(TSet)}) then                 error "The temperatures must be unique.  Two values are the same."                end if;   end proc;   TLIH:=proc() description "update LIH azeotrope";              if nargs > 0 then                 if type(args[1],`numeric`) then TLIHval:=args[1]                 else error "Expected a numeric value,%1",args[1] end if              else TLIHval:=NULL end if;              TSet:=[TLval,TIval,THval,TLIval,TIHval,TLHval,TLIHval];              if TLIHval <> NULL then printf("TLIH = %g",TLIHval) end if;              if nops(TSet)>nops({op(TSet)}) then                 error "The temperatures must be unique.  Two values are the same."                end if;   end proc;   checkT:=proc()           description "checks for valid order of TL, TI and TH";           if TLval <> NULL and TIval <> NULL and TLval > TIval then    printf("Bad temperature values; TL = %g and TI = %g", TLval,TIval) end if;           if TLval <> NULL and THval <> NULL and TLval > THval then printf("Bad temperature values: TL = %g and TH = %g", TLval, THval) end if;           if TIval <> NULL and THval <> NULL and TIval > THval then printf("Bad temperature values: TI = %g and TH = %g", TIval, THval) end if;   end proc;   MapSelectionNumber:=proc(Number::posint)        description "select the desired map given multiple feasible maps";        if nargs = 0 then MapSelection:=1        elif Number > 3 then           error "Map selection must be between 1 and 3"        else  MapSelection:=Number        end if   end proc;           DBM:=proc(Tseq)        description "find the DBM(s) given a sequence";        local soln, pick;        if whattype(TSequences[Tseq])=`indexed` then           error "The temperature sequence is invalid: %1.", Tseq end if;        soln:=[TSequences[Tseq]];        pick:=MapSelection;        if pick>nops(soln) then pick:=nops(soln) end if;        if nops(soln)>1 then           printf("The valid maps are: %s %s %s \n",op(soln));           printf("Selected map: %s",soln[pick]);        end if;        soln[pick];   end proc;      printData:=proc()              description "print the position temperature data";              if TLval <> NULL then printf("TL = %g \n", TLval) end if;              if TIval <> NULL then printf("TI = %g \n", TIval) end if;              if THval <> NULL then printf("TH = %g \n", THval) end if;              if TLIval <> NULL then printf("TLI = %g \n ",TLIval) end if;              if TIHval <> NULL then printf("TIH = %g \n ",TIHval) end if;              if TLHval <> NULL then printf("TLH = %g \n ",TLHval) end if;              if TLIHval <> NULL then printf("TLIH = %g \n ",TLIHval) end if;   end proc;   TSequence:=proc()     description "determine the temperature sequence from the temp order";     local tset,Fsort;     if TLval = NULL or TIval = NULL or THval = NULL then       printf("Missing key data; enter the component temperatures.");       return FAIL;     end if;     tset:= [TLval,1],[TIval,3],[THval,5];     if TLIval<>NULL then tset:=tset,[TLIval,2] end if;     if TIHval<>NULL then tset:=tset,[TIHval,4] end if;     if TLHval<>NULL then tset:=tset,[TLHval,6] end if;     if TLIHval<>NULL then tset:=tset,[TLIHval,7] end if;     tset:=[tset];     Fsort:=(x,y)->if op(1,x)<op(1,y) then true else false end if;     tset:=sort(tset,Fsort);     op(sscanf(cat(seq(convert(tset[i][2],`string`),i=1..nops(tset))), "%d"));   end proc;   # plot the ternary map   PlotDBM:=proc()     description "display the ternary map with all its data";     try       plotNodeSaddle(); # calculate the nodes and saddles       plotBoundaries(); # plot the boundary lines       plots[display](p1,p2,p3,p4); # display all of the plots together     end try;   end proc;   AzLI:=proc(xL::numeric)       description "enter the composition of the LI azeotrope";       if nargs=0 then AzLIval:=NULL       else         if xL > 0 and xL < 1 then            AzLIval:=xL;            printf("AzLI composition = %g in the low boiler",AzLIval);         else error "enter the low boiler composition as a value between 0 and 1."  end if;       end if;     end proc;     AzIH:=proc(xI::numeric)       description "enter the composition of the IH azeotrope";       if nargs=0 then AzIHval:=NULL       else         if xI > 0 and xI < 1 then            AzIHval:=xI;            printf("AzIH composition = %g in the intermediate boiler", AzIHval);         else error "enter the intermediate boiler composition as a value between 0 and 1."  end if;       end if;     end proc;   AzLH:=proc(xL::numeric)       description "enter the composition of the LH azeotrope";       if nargs=0 then AzLHval:=NULL       else         if xL > 0 and xL < 1 then            AzLHval:=xL;            printf("AzLH composition = %g in the low boiler", AzLHval);         else error "enter the low boiler composition as a value between 0 and 1."  end if;       end if;     end proc;   AzLIH:=proc(xL::numeric,xI::numeric)       description "enter the composition values for the ternary azeotrope";       if nargs=0 then AzLIHvalL:=NULL; AzLIHvalI:=NULL;       elif nargs=1 then error "enter two compositions for the low and intermediate boiling compounds"       else         if xL > 0 and xL < 1 and xI > 0 and xI < 1 and xL+xI<1 then            AzLIHvalL:=xL;            AzLIHvalI:=xI;            printf("AzLIH composition = %g in the low boiler and %g in the intermediate boiler", AzLIHvalL, AzLIHvalI);         else error "enter the compositions as values between 0 and 1.  Their sum must be less than 1."                  end if;       end if;     end proc;   MoveToTriangular:=(x,y,xm,ym)->op(evalf([x+y*cos(Pi/3)+xm, y*cos(Pi/6)+ym]));   plotNodeSaddle:=proc() local tseq, ns1, ns2, ns3, ns4, ns5, ns6, ns7, o1, o2, o3, o4, o5, o6, o7, i, j, out1, out2, Lpoint, Ipoint, Hpoint, LIpoint, IHpoint, LHpoint, LIHpoint, map, x, y;    description "create plot structure for showing the nodes and saddles on the DBM";    p3 := NULL;    tseq := convert(T:-TSequence(),`string`);  # use proc for sequence    for i from 1 to 7 do        o||i:= searchtext(convert(i,`string`),tseq)    end do;    # check position 1 for type (unstable node, stable node or saddle)    out1:=1; out2:=1; # start with assumption that T rises out of position 1    if o||2>0 and o||2<o||1 then out1:=0 end if;    if o||6>0 and o||6<o||1 then out2:=0 end if;    Lpoint:="UN"; # assume position 1 is an unstable node    if out1<>out2 then Lpoint:="S" end if; # saddle    if out1+out2=0 then Lpoint:="SN" end if; # stable node    # check position 3 for type (unstable node, stable node or saddle)    out1:=0; out2:=1;    if o||2>0 and o||3<o||2 then out1:=1 end if;    if o||4>0 and o||4<o||3 then out2:=0 end if;    Ipoint:="UN"; # assume position 3 is an unstable node    if out1<>out2 then Ipoint:="S" end if; # saddle    if out1+out2=0 then Ipoint:="SN" end if; # stable node    # check position 5 for type (unstable node, stable node or saddle)    out1:=0; out2:=0;    if o||4>0 and o||5<o||4 then out1:=1 end if;    if o||6>0 and o||5<o||6 then out2:=1 end if;    Hpoint:="UN"; # assume position 3 is an unstable node    if out1<>out2 then Hpoint:="S" end if; # saddle    if out1+out2=0 then Hpoint:="SN" end if; # stable node    map:=T:-DBM(Sequence);    # check position 2 for type (NULL, unstable node, stable node or saddle)    if map[1] = "0" then LIpoint = NULL    elif map[1] = "1" then LIpoint:="UN"    elif map[1] = "2" then LIpoint:="S"    elif map[1] = "3" then LIpoint:="SN"    else  LIpoint:="S" end if;    # check position 4 for type (NULL, unstable node, stable node or saddle)    if map[2] = "0" then IHpoint:=NULL    elif map[2] = "1" then IHpoint:="UN"    elif map[2] = "2" then IHpoint:="S"    elif map[2] = "3" then IHpoint:="SN"    else IHpoint = "S" end if;    # check position 6 for type (NULL, unstable node, stable node or saddle)    if map[3] = "0" then LHpoint:=NULL    elif map[3] = "1" then LHpoint:="UN"    elif map[3] = "2" then LHpoint:="S"    elif map[3] = "3" then LHpoint:="SN"    else LHpoint:="S" end if;    # check position 7 for type (NULL, unstable node, stable node or saddle)    if length(map) = 3 then LIHpoint:=NULL    elif map[5] = "m" then LIHpoint:="UN"    elif map[5] = "M" then LIHpoint:="SN"    else LIHpoint:="S" end if;    # check that the azeotrope data are available    if not(LIpoint=NULL) and AzLI = NULL then error("Enter AzLI") fi;    if not(IHpoint=NULL) and AzIH = NULL then error("Enter AzIH") fi;    if not(LHpoint=NULL) and AzLH = NULL then error("Enter AzLH") fi;    if not(LIHpoint=NULL) and (AzLIHvalL = NULL or AzLIHvalI = NULL) then error("Enter AzLIH") end if;    # plot position 1, L    x:=0; y:=1;  # in the plot, x = heavy boiler & y = low boiler    if Lpoint = "UN" then ns1:=rectangle([MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,0.01)], color=`blue`, scaling=`CONSTRAINED`)    elif Lpoint = "SN" then ns1:=rectangle([MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,0.01)], color=`red`, scaling=`CONSTRAINED`)    else ns1:=polygon([[MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,-0.01)], [MoveToTriangular(x,y,0,0.01)]], color=`turquoise`, scaling=`CONSTRAINED`)    end if;    # plot position 2, LI    x:=0; y:=AzLIval;  # in the plot, x = heavy boiler & y = low boiler    if LIpoint = "UN" then ns2:=rectangle([MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,0.01)], color=`blue`, scaling=`CONSTRAINED`)    elif LIpoint = "SN" then ns2:=rectangle([MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,0.01)], color=`red`, scaling=`CONSTRAINED`)    else ns2:=polygon([[MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,-0.01)], [MoveToTriangular(x,y,0,0.01)]], color=`turquoise`, scaling=`CONSTRAINED`)    end if;    # plot position 3, I    x:=0; y:=0;  # in the plot, x = heavy boiler & y = low boiler    if Ipoint = "UN" then ns3:=rectangle([MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,0.01)], color=`blue`, scaling=`CONSTRAINED`)    elif Ipoint = "SN" then ns3:=rectangle([MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,0.01)], color=`red`, scaling=`CONSTRAINED`)    else ns3:=polygon([[MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,-0.01)], [MoveToTriangular(x,y,0,0.01)]], color=`turquoise`, scaling=`CONSTRAINED`)    end if;    # plot position 4, IH    x:=1-AzIHval; y:=0;  # in the plot, x = heavy boiler & y = low boiler    if IHpoint = "UN" then ns4:=rectangle([MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,0.01)], color=`blue`, scaling=`CONSTRAINED`)    elif IHpoint = "SN" then ns4:=rectangle([MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,0.01)], color=`red`, scaling=`CONSTRAINED`)    else ns4:=polygon([[MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,-0.01)], [MoveToTriangular(x,y,0,0.01)]], color=`turquoise`, scaling=`CONSTRAINED`)    end if;    # plot position 5, H    x:=1; y:=0;  # in the plot, x = heavy boiler & y = low boiler    if Hpoint = "UN" then ns5:=rectangle([MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,0.01)], color=`blue`, scaling=`CONSTRAINED`)    elif Hpoint = "SN" then ns5:=rectangle([MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,0.01)], color=`red`, scaling=`CONSTRAINED`)    else ns5:=polygon([[MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,-0.01)], [MoveToTriangular(x,y,0,0.01)]], color=`turquoise`, scaling=`CONSTRAINED`)    end if;    # plot position 6, LH    x:=1-AzLHval; y:=AzLHval;  # in the plot, x = heavy boiler & y = low boiler    if LHpoint = "UN" then ns6:=rectangle([MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,0.01)], color=`blue`, scaling=`CONSTRAINED`)    elif LHpoint = "SN" then ns6:=rectangle([MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,0.01)], color=`red`, scaling=`CONSTRAINED`)    else ns6:=polygon([[MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,-0.01)], [MoveToTriangular(x,y,0,0.01)]], color=`turquoise`, scaling=`CONSTRAINED`)    end if;    # plot position 7, LIH    x:=1-AzLIHvalL-AzLIHvalI; y:=AzLIHvalL;  # in the plot, x = heavy boiler & y = low boiler    if LIHpoint = "UN" then ns7:=rectangle([MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,0.01)], color=`blue`, scaling=`CONSTRAINED`)    elif LIHpoint = "SN" then ns7:=rectangle([MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,0.01)], color=`red`, scaling=`CONSTRAINED`)    else ns7:=polygon([[MoveToTriangular(x,y,-0.01,-0.01)], [MoveToTriangular(x,y,0.01,-0.01)], [MoveToTriangular(x,y,0,0.01)]], color=`turquoise`, scaling=`CONSTRAINED`)    end if;    # prepare the combined plot of nodes and saddles    p3:=plots[display](ns1,ns2,ns3,ns4,ns5,ns6,ns7);   end proc;   plotBoundaries:=proc()     description "draw straight lines for distillation boundaries";     local tseq, dbm, boundaries, getpoint, i, position1, position2, line;     p4:=NULL;     line:=NULL;     getpoint:=proc(position::`string`)       description "create x-y point for the position";       if position = "1" then          [0,1]       elif position = "2" then          [0,AzLIval]       elif position = "3" then          [0,0]       elif position = "4" then          [1-AzIHval,0]       elif position = "5" then          [1,0]       elif position = "6" then          [1-AzLHval,AzLHval]       else          [1-AzLIHvalL-AzLIHvalI,AzLIHvalL]       end if;     end proc;     try       tseq:=TSequence();       dbm:=DBM(tseq);       boundaries:=eval(Boundaries[dbm]);     end try;     if nops(boundaries) = 0 then 0     else       for i from 1 to nops(boundaries) do         position1,position2:=op(sscanf(convert(boundaries[i], `string`),"%c%c"));         line:=line,[getpoint(position1),getpoint(position2)];       end do;       p4:=plot({line}, scaling=`CONSTRAINED`, color=`black`, axes=`NONE`, coords=`triangular`,thickness=2);       nops(boundaries);     end if;   end proc;   AddPoint:=proc()     description "add a point to the points set";     local makept, xL, xI, xH, ptname, firstpoint, secondpoint, MixType, MixSum, amt, slope1, slope2, b1, b2;     makept:=proc(f::`name`,nam,typ,xLi,xIi)         description "make a composition point record";         f:=module() export name, type, xL, xI; name:=nam; type:=typ; xL:=xLi; xI:=xIi; end module;         f;     end proc;     if nargs=0 then error "arguments are required"     elif nargs=1 then       if type(args[1],`string`) then          print("deleting the composition point");          Raw_Materials:=Raw_Materials minus {convert(args[1],`name`)};          Products:=Products minus {convert(args[1],`name`)};          points:=points minus {convert(args[1],`name`)};       else  error "wrong type of argument 1; string required; received %1", args[1]       end if;     elif nargs=2 then       if type(args[1],`string`) then          if member(args[2], {"Feed","Intermediate","Product"}) then             xL:=0;             xI:=0.5;             ptname:=convert(args[1],`name`);             if args[2] = "Feed" then               Raw_Materials:=Raw_Materials union {ptname}             elif args[2] = "Product" then               Products:=Products union {ptname}             end if;             points:=points minus {ptname};             points:=points union {makept(ptname, args[1], args[2] ,xL, xI)};          else            error "error in argument 2; it must be \"Feed\", \"Intermediate\", or \"Product\"; received %1", args[2]          end if       else error "wrong type of argument 1; string required; received %1", args[1]       end if;     elif nargs=4 then       if not type(args[1],`string`) then          error "a string is required for argument 1 to enter point name"  end if;       if not member(args[2], {"Feed","Intermediate","Product"}) then error "the second argument must be \"Feed\", \"Intermediate\" or \"Product\" to tell the type of point"  end if;       if not (type(args[3],`numeric`) or type(args[3],`+`)) then error "the third argument must be numeric for xL or an expression like A+B for mixing existing points" end if;       if not type(args[4],`numeric`) then error "the fourth argument must be numeric" end if;       if type(args[3],`numeric`) then             xL:=args[3];             xI:=args[4];             ptname:=convert(args[1],`name`);             if args[2] = "Feed" then               Raw_Materials:=Raw_Materials union {ptname}             elif args[2] = "Product" then               Products:=Products union {ptname}             end if;             points:=points minus {ptname};             points:=points union {makept(ptname, args[1],  args[2], xL, xI)};       else         firstpoint:=convert(op(1,args[3]),`name`);         secondpoint:=convert(op(2,args[3]),`name`);         MixType:=1;         if nops(op(2,args[3]))=2 then            secondpoint:=op(2,op(2,args[3]));            secondpoint:=convert(secondpoint,`name`);            MixType:=-1;         end if;         if not member(firstpoint, points) then               error "%1 of %2 is not a defined point", op(1,args[3]), args[3]         end if;         if not member(secondpoint, points) then               error "%1 of %2 is not a defined point", op(2,args[3]), args[3]         end if;         if MixType=1 then            xL:=firstpoint:-xL+secondpoint:-xL*args[4];            xI:=firstpoint:-xI+secondpoint:-xI*args[4];            xH:=1. - firstpoint:-xL - firstpoint:-xI + (1. - secondpoint:-xL - secondpoint:-xI) * args[4];            MixSum:=xL+xI+xH;            xL:=xL/MixSum;            xI:=xI/MixSum;         else            amt:=args[4];            xL:=firstpoint:-xL-secondpoint:-xL*amt;            if xL<0 then               amt:=firstpoint:-xL/secondpoint:-xL;               xL:=0;            end if;            xI:=firstpoint:-xI-secondpoint:-xI*amt;            if xI<0 then               amt:=firstpoint:-xI/secondpoint:-xI;               xI:=0;            end if;            xH:=1. - firstpoint:-xL - firstpoint:-xI - (1. - secondpoint:-xL - secondpoint:-xI) * amt;            if xH<0 then               amt:= (1. - firstpoint:-xL - firstpoint:-xI) / (1. - secondpoint:-xL - secondpoint:-xI);               xH:=0;            end if;            xL:=firstpoint:-xL-secondpoint:-xL*amt;            xI:=firstpoint:-xI-secondpoint:-xI*amt;            MixSum:=xL+xI+xH;            xL:=xL/MixSum;            xI:=xI/MixSum;         end if;         ptname:=convert(args[1],`name`);         if args[2] = "Feed" then               Raw_Materials:=Raw_Materials union {ptname}         elif args[2] = "Product" then               Products:=Products union {ptname}         end if;         points:=points minus {ptname};         points:=points union {makept(ptname, args[1],  args[2], xL, xI)};       end if     elif nargs=3 then       if not type(args[1],`string`) then  error "a string is required for argument 1 to enter point name"  end if;       if not member(args[2], {"Feed","Intermediate","Product"}) then error "the second argument must be \"Feed\", \"Intermediate\" or \"Product\" to tell the type of point"  end if;       if not type(args[3],`=`) then error "the third argument must be a formula of two points equal to two other points (A+B=C+D)" end if;       if not (nops(op(1,args[3]))=2 and nops(op(2,args[3]))=2) then error "the formula for argument 3 must be in the form of two points = two points" end if;       firstpoint:=convert(op(1,op(1,args[3])),`name`);       secondpoint:=convert(op(2,op(1,args[3])),`name`);       if not member(firstpoint, points) then   error "%1 of %2 is not a defined point", firstpoint, args[3]  end if;       if not member(secondpoint, points) then  error "%1 of %2 is not a defined point", secondpoint, args[3]  end if;       firstpoint:=convert(op(1,op(2,args[3])),`name`);       secondpoint:=convert(op(2,op(2,args[3])),`name`);       if not member(firstpoint, points) then   error "%1 of %2 is not a defined point", firstpoint, args[3]  end if;       if not member(secondpoint, points) then  error "%1 of %2 is not a defined point", secondpoint, args[3] end if;       slope1:=(secondpoint:-xI - firstpoint:-xI) / (secondpoint:-xL - firstpoint:-xL);       b1:=firstpoint:-xI - slope1*firstpoint:-xL;       firstpoint:=convert(op(1,op(1,args[3])),`name`);       secondpoint:=convert(op(2,op(1,args[3])),`name`);       slope2:=(secondpoint:-xI - firstpoint:-xI) / (secondpoint:-xL - firstpoint:-xL);       b2:=firstpoint:-xI - slope2*firstpoint:-xL;       xL:=(b2 - b1) / (slope1 - slope2);       xI:=slope1*xL + b1;       ptname:=convert(args[1],`name`);       if args[2] = "Feed" then            Raw_Materials:=Raw_Materials union {ptname}       elif args[2] = "Product" then            Products:=Products union {ptname}       end if;       points:=points minus {ptname};       points:=points union {makept(ptname, args[1],  args[2], xL, xI)};     end if   end proc;   PrintPoints:=proc()     description "prints data on all defined points";     local i;     print("pt name, pt type, xL, xI");     for i in points do       printf("s s 8.5f 8.5f",i:-name, i:-type, evalf(i:-xL), evalf(i:-xI))     end do;   end proc;   AddUnit:=proc(ID::`string`,Unit::`equation`)     description "add a unit to the toolbox";     local i, input, output, unitname, makeunit;     if nargs=1 then       print("deleting the unit");       Units:=Units minus { convert(ID,`name`) };       return Units;       end if;        makeunit:=proc(g::`name`,ID,ins,outs)         description "make a unit op record";         g:=module() export name, inputs, outputs; name:=ID; inputs:=ins; outputs:=outs; end module;         g;     end proc;     if nops(lhs(Unit)) = 1 then        input:=lhs(Unit);        if not member(input,points) then           error "%1 is not a defined point", input        end if;     else        input:={op(lhs(Unit))};        input:=seq(convert(input[i],`name`),i=1..nops(input));        for i in {op(lhs(Unit))} do            if not member(i, points) then               error "%1 is not a defined point", i            end if        end do;     end if;     if nops(rhs(Unit)) = 1 then        output:=rhs(Unit);        if not member(output,points) then           error "%1 is not a defined point", output        end if;     else        output:={op(rhs(Unit))};        output:=seq(convert(output[i],`name`),i=1..nops(output));        for i in {op(rhs(Unit))} do            if not member(i, points) then               error "%1 is not a defined point", i            end if        end do;     end if;     unitname:=convert(ID,`name`);     Units:=Units minus {unitname};     Units:=Units union {makeunit(unitname,ID,{input},{output})};     Units;   end proc;       PrintPoints:=proc()     description "prints data on all defined points";     local i;     print("pt name, pt type, xL, xI");     for i in points do       print(i:-name, i:-type, evalf(i:-xL), i:-xI)     end do;   end proc;       Psi_minus:=proc(OpSet)            description "Psi- function of combinatorial analysis";            local i, solution;            solution:={};            for i from 1 to nops(OpSet) do              solution:=solution union op(1,op(i,OpSet))            end do;            solution;   end proc:   Psi_plus:=proc(OpSet)           description "Psi+ function of combinatorial analysis method";           local i, solution;           solution:={};           for i from 1 to nops(OpSet) do             solution:=solution union op(2,op(i,OpSet))           end do;           solution;   end proc:     phi_minus:=proc(mSet,OpSet)            local i, solution, item;            description "phi- function of combinatorial analysis method";            solution:={};            for i from 1 to nops(OpSet) do                item:={op(i,OpSet)};                if (op(2,op(1,item)) intersect mSet) <> {} then                   solution:=solution union item  end if;            end do;;            solution;   end proc:   phi_plus:=proc(mSet,OpSet)           description "phi+ function of combinatorial method";           local i, solution, item;           solution:={};           for i from 1 to nops(OpSet) do             item:={op(i,OpSet)};             if (op(1,op(1,item)) intersect mSet) <> {} then                solution:=solution union item  end if;           end do;;           solution;   end proc:   MSG:=proc()      description "MSG algorithm for flowsheet maximal structure";      local i, item, opset, Mat, o, r, p, m, ox;      Op:={seq([Units[i]:-inputs,Units[i]:-outputs],i=1..nops(Units))};      opset:=Op minus phi_minus(Raw_Materials,Op);      Mat:=Psi_minus(opset) union Psi_plus(opset);      r:=Psi_minus(opset) minus (Psi_plus(opset) union Raw_Materials);      if nops(r) > 0 then         while r <> {}  do               item:={op(1,r)};               Mat:=Mat minus item;               o:=phi_plus(item,opset);               opset:=opset minus o;               r:=(r union (Psi_plus(o) minus Psi_plus(opset))) minus item;         end do;      end if;      if (Products union Mat) <> Mat then         error "No maximal structure exists!"      end if;      p:=Products;      m:={};      o:={};      while p <> {}  do         item:={op(1,p)};         m:=m union item;         ox:=phi_minus(item,opset);         o:=o union ox;         p:=(p union Psi_minus(ox)) minus (Raw_Materials union m);      end do;      m:=Psi_minus(o) union Psi_plus(o);      maximal_structure:=[m,o];      maximal_structure;   end proc:   graph_network:=proc(Structure)    description "convert a Pgraph structure into a Maple network";    local i, j, rm, pr, mid, m, mat_order, listing, opset, Gr, mat_ins, mat_outs, el, unitnum, inset, outset, inset2, outset2, unit_order, unit_order2, unitname, location;    m:=op(1,Structure);    opset:=op(2,Structure);    # check for material inlets and outlets of the structure    mat_ins:=Psi_minus(opset);    mat_outs:=Psi_plus(opset);    # find all raw materials (they are only inlets to unit ops)    rm:={};    # find all products (they are only outlets of unit ops)    pr:={};    for el in m  do        if member(el,mat_ins) and member(el,mat_outs) = false then           rm:=rm union {el} end if;        if member(el,mat_ins) = false and member(el,mat_outs) then           pr:=pr union {el} end if;    end do;    mid:=m minus rm minus pr;    mat_order:=[op(pr)],[op(mid)],[op(rm)];    # get labels from Units information    j:=0;    for el in opset do        j:=j+1;        inset:=op(1,el);        outset:=op(2,el);        for i from 1 to nops(Units) do            inset2:=Units[i]:-inputs;            outset2:=Units[i]:-outputs;            if inset = inset2 and outset = outset2 then break fi;        end do;        unitname:=Units[i]:-name;        unitnum[j]:=unitname;    end do;    unit_order:=[seq(cat(`:`,unitnum[j]), j=1..nops(opset))];    for i from 1 to nops(unit_order) do        if searchtext("_",convert(unit_order[i],`string`)) > 0 then           unit_order[i]:= convert( substring( convert( unit_order[i],  `string`), 1 .. searchtext("_", convert(unit_order[i], `string`))-1), `name`);        end if;    end do;    unit_order2:=convert(convert(unit_order,`set`),`list`);    listing:=mat_order,unit_order2;    new(Gr);    addvertex(op(1,Structure),Gr);    addvertex(unit_order,Gr);    for i from 1 to nops(opset) do        connect(op(1,op(i,opset)), [op(i,unit_order)], directed, Gr);        connect([op(i,unit_order)], op(2,op(i,opset)), directed, Gr);    end do;    Linear(listing),Gr;   end proc:   SSG:=proc(by_prod)      description "SSG algorithm to find all feasible flowsheets";    local mats, e;      by_prods:=by_prod;      mats:=Psi_minus(Op) union Psi_plus(Op);    ox_vals:=table();      for e in mats do           ox_vals[{e}]:=phi_minus({e},Op);      end do;    SSGmethod(Products,{},{});      NumberOfMaps:=nops(decision_map);    print("# of decision maps found = "||NumberOfMaps);      decision_map;   end proc:   SSGmethod:=proc(p,m,delta_m)    description "supporting function for SSG method";    local i, j, k, item_x, set_c, testing, item_c, item_y, ox, oy, delta_y, delta_m2, p2, m2, op_y, mat_outs, mat_ins, opmat, pr, el;    if nops(p) = 0 then      if by_prods = true then        # by-products are allowed; accept the solution        decision_map:=decision_map union {delta_m}      elif by_prods = false then        # accept a solution if no by-products are present        # check for material inlets and outlets of the decision map        opmat:={seq(op(op(2,op(k,delta_m))),k=1..nops(delta_m))};        mat_ins:=Psi_minus(opmat);        mat_outs:=Psi_plus(opmat);        # subtract the raw materials and products        mat_ins:=mat_ins minus Raw_Materials minus Products;        mat_outs:=mat_outs minus Raw_Materials minus Products;        # for a valid mapping, the intermediate materials must be both        # inlet and outlet materials        if ( (mat_ins minus mat_outs) = {} and    (mat_outs minus mat_ins) = {} ) then           decision_map:=decision_map union {delta_m};        end if;      else        # accept the solution if the by-products are within the input set        # of allowed by-products.        # check for material inlets and outlets of the structure        opmat:={seq(op(op(2,op(k,delta_m))),k=1..nops(delta_m))};        mat_ins:=Psi_minus(opmat);        mat_outs:=Psi_plus(opmat);        # find all products (they are only outlets of unit ops)        pr:={};        for el in (mat_ins union mat_outs) do           if member(el,mat_ins) = false and member(el,mat_outs) then              pr:=pr union {el} end if;        end do;        # test that the by-products are contained in the by_prods set.        # accept the solution if true.        if ((pr minus Products) minus by_prods) = {} then           decision_map:=decision_map union {delta_m} end if;      end if;      if (nops(decision_map) mod 100) = 0 and          nops(decision_map) > 0              then          print(nops(decision_map), ` decision maps found so far`)      end if;      return 0;    end if;    item_x:={op(1,p)};    ox:=ox_vals[item_x];    set_c:=choose(ox) minus {{}};    for i from 1 to nops(set_c) do        item_c:=op(i,set_c);        testing:=0;        for j from 1 to nops(m) do            item_y:={op(j,m)};            oy:=ox_vals[item_y];            delta_y:={};            if nops(delta_m) > 0 then               op_y:={seq(op(op(2,op(k,delta_m))),k=1..nops(delta_m))};               delta_y:=phi_minus(item_y,op_y);            end if;            if (item_c intersect (oy minus delta_y)) <> {}  or                  ((ox minus item_c) intersect delta_y) <> {}  then                 testing:=1;                 break;            end if;        end do;        if testing = 1  then  next  end if;        delta_m2:=delta_m union {[item_x,item_c]};        p2:=p;        for j from 1 to nops(item_c) do            p2:=p2 union op(1,op(j,item_c))        end do;        p2:=p2 minus (Raw_Materials union m union item_x);        m2:=m union item_x;        SSGmethod(p2, m2, delta_m2);    end do;   end proc:     flowsheet:=proc(item::integer)     description "draw a simple network graph of a selected flowsheet";     if maximal_structure = {} then        error "the maximal structure from procedure MSG is not available; run MSG."     end if;     if decision_map = {} then        error "the decision map from procedure SSG is not available; run SSG."     end if;     if nops(decision_map) < item then        error "there are only %1 flowsheet(s) in the decision map and you entered %2.", nops(decision_map), item     end if;     if item < 0 then        error "enter a positive integer"     end if;     if item = 0 then        draw(graph_network(maximal_structure))     else        draw(graph_network(Pgraph(decision_map[item])))     end if;   end proc;     Pgraph:=proc(delta_m)     description "convert a given decision map into a P graph";     local k, m, o;     o:={seq(op(op(2,op(k,delta_m))),k=1..nops(delta_m))};     m:=Psi_minus(o) union Psi_plus(o);     [m,o];   end proc;   PlotPoints:=proc(item::integer)     description "plot points from a flowsheet on the DBM";     local m, o, i, diskcolor, eqn, eqn2, j, point1, point2, titleinfo;     eqn2:="Units: ";     if item < 0 then error "the first argument is %1; it must be 0 or greater", item end if;     if maximal_structure = {} then        error "the maximal structure from procedure MSG is not available; run MSG."     end if;     if decision_map = {} then        error "the decision map from procedure SSG is not available; run SSG."     end if;     if nops(decision_map) < item then        error "there are only %1 flowsheet(s) in the decision map and you entered %2.", nops(decision_map), item     end if;     if item = 0 then       m,o:=op(maximal_structure);       titleinfo:="Flowsheet for Maximal Structure";     else       m,o:=op(Pgraph(decision_map[item]));       titleinfo:="Flowsheet #"||item;     end if;     for i from 1 to nops(m) do         if m[i]:-type = "Feed" then diskcolor:=`green`         elif m[i]:-type = "Product" then diskcolor:=`navy`         else diskcolor:=`gold`  end if;         pt||i:=disk([MoveToTriangular(1.- m[i]:-xL - m[i]:-xI, m[i]:-xL, 0,0)], 0.018, color=diskcolor, scaling=`CONSTRAINED`)     end do;     for i from 1 to nops(m) do         if m[i]:-type = "Feed" then diskcolor:=`green`         elif m[i]:-type = "Product" then diskcolor:=`navy`         else diskcolor:=`gold`  end if;         t||i:=plots[textplot]( [ MoveToTriangular(1. - m[i]:-xL - m[i]:-xI, m[i]:-xL, 0.035, 0), m[i]:-name], color=diskcolor);     end do;     p5:=plots[display](seq(pt||i,i=1..nops(m)), seq(t||i, i=1..nops(m)) );     j:= 0;     for i from 1 to nops(o) do         if nops(o[i][1]) = 2 then           j:=j+1;           point1:=o[i][1][1];           point2:=o[i][1][2];           l||j:=line([MoveToTriangular(1. - point1:-xL - point1:-xI, point1:-xL,0,0)],[MoveToTriangular(1. - point2:-xL - point2:-xI, point2:-xL, 0,0)],color=`red`,linestyle=3,thickness=1);         end if;         if nops(o[i][2]) = 2 then           j:=j+1;           point1:=o[i][2][1];           point2:=o[i][2][2];           l||j:=line([MoveToTriangular(1. - point1:-xL - point1:-xI, point1:-xL,0,0)],[MoveToTriangular(1. - point2:-xL - point2:-xI, point2:-xL, 0,0)],color=`red`,linestyle=1,thickness=1);         end if;           end do;     p6:=plots[display](seq(l||i,i=1..j));     for i from 1 to nops(o) do       eqn:=convert(cat(convert(convert( f(op(o[i][1])),`+`), `string`), "=", convert(convert( f(op(o[i][2])), `+`), `string`) ), `name` );       eqn2:=cat(eqn2,eqn,"; ");     end do;     MapUnits:=eqn2;     plots[display](p6, p1, p2, p3, p4, p5,title=titleinfo);     end proc; end module; # updating the draw package to show different colors and line styles on the graph `draw/Linear`:= proc (partitions::specfunc(list,Linear), G::GRAPH, origin::list(numeric), xrng::name, yrng::name)   local e, max_x, middle, n, nx, part, shift, t1, text, pos, x, y, t, vset, lines, points, i, j, v, center, a, b, c, d, a1, a2, a3, angle, dist; # add new local vars   option `Copyright (c) 1992 by the University of Waterloo. All rights reserved.`;   center := origin[1 .. 2];   vset := {};   n := 0;   part := table();   shift := table();   for t in partitions while true do     n := n+1;     if not type(t,list) then error "partition should be a list" end if;     shift[n] := 0;     t1 := select(hastype,t,('identical')('offset') = numeric);     if t1 <> [] then       shift[n] := subs({op(t1)},'offset');       t := select(proc (x) not hastype(x,`=`) end proc,t)     else if `intersect`({op(t)},vset) <> {} then       error "intersecting partitions involving %1", t end if     end if;     if not type(t,('list')(('VERTEX')(G))) then         error "not a list of vertices %1", t     end if;     part[n] := t;     vset := `union`(vset,{op(t)})   end do;   userinfo(5,'networks',`allocated vertices:`,print(vset));   if `minus`(networks['vertices'](G),vset) <> {} then     n := n+1;     t := sort([op(`minus`(networks['vertices'](G),vset))]);     part[n] := t;     shift[n] := 0;     userinfo(1,'networks',`adding a new partition`,t)   end if;   max_x := origin[1];   for j to n do     if max_x < nops(part[j]) then max_x := nops(part[j]) end if   end do;   middle := 1/2*max_x-1/2+origin[1];   points := table();   text := table();   pos := table();   x := origin[1];   for j to n do     vset := part[j];     nx := nops(vset);     x := x+1;     # update it to shift the alignment of columns     y := middle-1/2*nx+1/2+shift[j]+origin[2]+0.2*((j mod 3)-1);     for i to nops(vset) do       v := vset[i];       pos[v] := [x, y];       points[v] := POINTS(pos[v]);       if j = n then text[v] := ('TEXT')([pos[v][1]+.8e-1,  pos[v][2]], convert(v,string))       else text[v] := ('TEXT')([pos[v][1]-.1, pos[v][2]], convert(v,string))       end if;       y := y+1     end do   end do;   lines := table();   for e in edges(G) while true do     x := networks['ends'](e,G);     if 1 < nops(x) then y := x[2]; x := x[1] else x := x[1] end if;     # update the line color and linestyle     if substring(convert(networks['head'](e,G),string),1..1) = ":"       then         lines[e] := CURVES([pos[x], pos[y]], COLOUR('RGB',1.0,0,0), LINESTYLE(0) )     else         lines[e] := CURVES([pos[x], pos[y]], COLOUR('RGB',0,0,1.0), LINESTYLE(3))     end if   end do;   t := map(op,{entries(points), entries(lines), entries(text)});   xrng := origin[1]+1/2 .. origin[1]+n+1/2;   yrng := origin[2]-1/2 .. origin[2]+max_x-1/2;   return t end proc:

Warning, these names have been redefined: dodecahedron, icosahedron, octahedron, tetrahedron

Warning, the protected name Chi has been redefined and unprotected

Warning, not an orthogonal coordinate system - no scale factors calculated.

Enter Ternary System Data

The system pressure is set and the ternary components are placed in boiling point temperature order:  TL, low boiler;  TI, intermediate boiler;  and TH, heavy boiler.  You enter the boiling point temperatures that apply for your system using the following procedures:  Enter nothing if the azeotrope does not exist. (i.e., T:-TLI(); )

>	T:-TL(45);  # low boiler T:-TI(55);  # intermediate boiler T:-TH(65);  # heavy boiler T:-TLI(41);   # L-I azeotrope (if it exists) T:-TIH(50);   # I-H azeotrope (if it exists) T:-TLH(42);   # L-H azeotrope (if it exists) T:-TLIH(25);  # ternary azeotrope (if it exists)

TL = 45       

TI = 55       

TH = 65       

TLI = 41       

TIH = 50       

TLH = 42       

TLIH = 25       

Determine the temperature sequence from the data using the T:-TSequence() procedure.  You can read more about temperature sequences from http://users.chartertn.net/lpartin.

>	Sequence:=T:-TSequence();

Select which map when there are multiple feasible maps (T:-DBM will tell you when there are multiple maps)

>	T:-MapSelectionNumber(1);

The software finds a simple form of the RCM called a distillation boundary map (DBM).  The DBM plot shows the boundaries that exist in the RCM.  If the DBM does not exist for your temperature sequence, then you must correct the temperature data from above.  Maybe you missed the existance of an azeotrope.

>	DBM:=T:-DBM(Sequence);

Enter the compositions for the azeotropes that exist:

>	T:-AzLI(0.88);  # enter L composition of LI azeotrope if it exists T:-AzIH(0.45);  # enter I composition of IH azeotrope if it exists T:-AzLH(0.8);  # enter L composition of LH azeotrope if it exists T:-AzLIH(0.45,0.25); # enter L and I compositions of the ternary azeotrope if it exists

AzLI composition = .88     in the low boiler

AzIH composition = .45     in the intermediate boiler

AzLH composition = .8      in the low boiler

AzLIH composition = .45     in the low boiler and .25     in the intermediate boiler

Plot the Ternary System

The DBM is plotted in triangular coordinates using the T:-PlotDBM() procedure.  Unstable nodes are marked as blue squares; they are the coolest temperature within their regions.  Stable nodes are marked as red squares; they are the highest temperatures within their regions.  Saddles are marked as turquoise triangles.  Distillation boundaries are shown as bold, straight lines.  

>	T:-PlotDBM();

Flowsheet Synthesis

The synthesis task is done by adding ternary composition points and then creating unit operations with the points as the input and output streams.  The software determines when a feasible maximal structure has been generated to make the defined products from the feeds (MSG algorithm).  It then finds all of the feasible flowsheets within the maximal structure (SSG algorithm).  You can plot a flowsheet as a simple network diagram or show the composition points on the ternary diagram.

Enter composition points and units

The AddPoint procedure is used to enter a composition point.  Call it once for each composition point.

  T:-AddPoint(pointName, pointType, xL, xI)
       pointName, a string giving the point name  (Do not use "e", "O", or "T" as point names.)

      pointType, a string giving the type of point: "Feed", "Intermediate" or "Product"

      xL, the point fractional composition in terms of the light (L) boiling component
      xI,  the point fractional composition in terms of the intermediate (I) boiling component

     (If only pointName is entered, then the point is deleted.  Be sure to also delete any unit that uses the point.)

Once you have two or more points defined, then you can use the existing points to define new points.  Two points can be mixed as follows:

  T:-AddPoint(pointName,pointType,mixDefinition,fractionSecondPoint)

      pointName, a string giving the point name  (Do not use "e", "O", or "T" as point names.)

      pointType, a string giving the type of point: "Feed", "Intermediate" or "Product"

      mixDefinition, definition of stream mixing as point1+point2.  Examples:  F+t, A+B, C+E

      fractionSecondPoint, a number giving the amount of point2 to add to one part of point1.

You can also create a new point by subtracting an existing point from another point as follows:

 T:-AddPoint(pointName,pointType,subtractDefinition,fractionSecondPoint)

      pointName, a string giving the point name  (Do not use "e", "O", or "T" as point names.)

      pointType, a string giving the type of point: "Feed", "Intermediate" or "Product"

      subtractDefinition, definition of stream subtraction as point1-point2.  Examples:  F-t, A-B, C-E

      fractionSecondPoint, a number giving the amount of point2 to subtract from one part of point1.

The final way to create is new point is to find the intersection of two lines.  Each line is defined by two points.

 T:-AddPoint(pointName,pointType,intersectionDefinition)

•    pointName, a string giving the point name  (Do not use "e", "O", or "T" as point names.)

      pointType, a string giving the type of point: "Feed", "Intermediate" or "Product"

      intersectionDefinition, four points in the form P1+P2=P3+P4.  The defined point is the intersection of the line from P1 to P2 with the line from P3 to P4.

>

T:-AddPoint("F","Feed",0.5,0.5): T:-AddPoint("W","Product",0,1.): T:-AddPoint("E","Product",1.,0): T:-AddPoint("t","Intermediate",0,0): T:-AddPoint("A","Intermediate",F+t,0.24): T:-AddPoint("B","Intermediate",0.375,0.50): T:-AddPoint("J","Intermediate",F+B,0.9): T:-AddPoint("C","Intermediate",A-B,0.77): T:-AddPoint("H","Intermediate",0.45,0.25): T:-AddPoint("K","Intermediate",J-W,0.3): T:-AddPoint("d","Intermediate",B-W,0.3): T:-AddPoint("M","Intermediate",C+K=E+H): T:-AddPoint("G","Intermediate",C-t,0.2): T:-AddPoint("o","Intermediate",d+G=E+H): T:-AddPoint("L","Intermediate",G+K=E+H);

>	# Printing the composition points T:-PrintPoints();

Unit operations are defined in terms of the composition points using AddUnit:

     T:-AddUnit(unitName, unitEquation);

            unitName, a string giving the name of the unit operation (i.e., "mix1")

            unitEquation, an equation of inputs = outputs in terms of the composition points  (i.e., A=B+C)
If you enter only the unitName, then the unit is deleted.

(Note:  You can combine several units into one unit when printing a network diagram of a flowsheet by giving the base name to the units followed by _a, _b, _c, etc.  For example, mix_a, mix_b and mix_c would plot as mix on the flowsheet plot.  This simplifies the flowsheet plot.)

>

T:-AddUnit("Mix1",F+t=A): T:-AddUnit("Decant1_a",A=B+C): T:-AddUnit("Distill1",B=d+W): T:-AddUnit("Distill2",C=G+t): T:-AddUnit("Mix2",G+d=o): T:-AddUnit("Distill3_a",o=H+E): T:-AddUnit("Decant1_b",H=B+C): T:-AddUnit("Mix3",F+B=J): T:-AddUnit("Distill4",J=K+W): T:-AddUnit("Mix4",G+K=L): T:-AddUnit("Distill3_b",L=H+E): T:-AddUnit("Mix5",C+K=M): T:-AddUnit("Distill3_c",M=H+E);

Solve for Feasible Flowsheets

First, the MSG algorithm finds the maximal superstructure if it exists.

>

T:-MSG();

If the maximal structure exists, then the SSG algorithm finds all feasible flowsheets.

>	T:-SSG(false):

View Flowsheet(s)

You can view a flowsheet as a simple network diagram via the flowsheet procedure.  It has a single argument that is the flowsheet number as found by the SSG procedure.  Enter flowsheet 0 for the maximal structure.

The information is placed in columns:

• products in first column
• intermediates in second column
• feeds in third column
• units in last column

Red lines connect the unit name to its feeds.  Dashed blue lines connect

>	T:-flowsheet(0);

You can view the flowsheet stream compositions on the ternary map using the T:-PlotPoints procedure.  It has a single argument which is the flowsheet number as found by the SSG procedure.  Enter flowsheet 0 for the maximal flowsheet structure.  

• Feed compositions are shown as green circles.
• Product compositions are navy circles.
• Intermediate compositions are gold circles.
• A dashed red line connects feed points when there are two feeds to a given unit.
• A solid red line connects product points when there are two products from a given unit.
• Stable nodes of the residue curve map are marked as red squares.
• Unstable nodes of the residue curve map are marked as blue squares.
• Saddle nodes of the residue curve map are marked as turquoise triangles.
• Distillation boundaries are shown as bold, straight lines between nodes.

>	T:-PlotPoints(1); T:-MapUnits;

Animate the plotting of all flowsheets

>	plots[display](seq(T:-PlotPoints(i),i=0..T:-NumberOfMaps),insequence=true);

Restart

The restart command resets all information.  You would then need to start at the beginning again by loading the T module code.

>

restart;

>