Thanks for the comments and feedback. Most of the code is in data preparation. The core of the methodology in matrix inversion and multiplication.
1) I will use the "TAB" instead of redefining it. 2) I am still in testing mode and so hard paths are OK at this time. If and when this is distributed to others, it should have a GUI, where user can selects the the required files. 3) jview is not needed for the script and safe to delete. 4) I inconsistently used global and local scope. I will make them local. 5) I checked the output against the published results from the reference I cited above the code and they agree. 6) This is a small test case. My intention is to apply to large to apply to large network where the matrix size can be 50000x50000 and this can be very sparse. PT On Fri, May 24, 2013 at 12:38 PM, Raul Miller <[email protected]> wrote: > --This is a lot of code and I do not have access to the reference > (that I know of), so my comments will be superficial. > > 1. TAB is already defined and has the same definition as your 'tab'. > > 2. Hardcoded paths are awkward. For something like this they are > probably fine. It's also nice that you included sample data to test > against. I had to fix some problems induced by email line wrap. > > 3. But I don't have jview, so I commented out the line that loads it > -- I don't know if that matters. > > 4. I got a domain error: > |domain error: script > | lineSusceptance =:(_1*+/lineSusceptance)setDiag lineSusceptance > > I think the issue here is that lineSusceptance was already defined > locally and here it's being defined globally. If you are running this > in a global context you will not have this problem, but loading it as > a script fails because of this distinction (you cannot update a global > name in the same context where that identical name exists as a local > name - this is easy to work around if you really want it to happen and > is almost always a coding mistake.) > > 5. I do not know how to verify that the results are correct, and am > generally not in a position to talk about correctness of the code. > > 6. Unless your data is much sparser than the example, I think that a > sparse representation will not help. > > Thanks, > > -- > Raul > > > On Fri, May 24, 2013 at 12:23 PM, P T <[email protected]> wrote: > > Here is the my first program in J! > > > > > > > > NB. Goal: To compute DC power flows in an electrical transmission network > > NB. given transmission line reactance and bus injections and withdrawals > > > > NB. Reference: R.D. CHRISTIE, B. WOLLENBERG, AND I. WANGENSTEEN > > NB. "Transmission Management in the Deregulated Environment" > > NB. Proceedings of IEEE, Vol. 88, N0. 2 Feb 2000 > > > > > > load 'files' > > load 'jview' > > > > tab=. 9{a. > > > > NB. read the bus and branch data from files > > > > NB. rawData=. <;._2 "1 (,.&tab 'm' fread 'E:\EMCAS Misc\J > > Loadflow\Bus11.txt') > > NB. bus=. _ ".&> rawData > > > > NB. rawData=. <;._2 "1 (,.&tab 'm' fread 'E:\EMCAS Misc\J > > Loadflow\Branch11.txt') > > NB. branch=. _ ".&> rawData > > > > NB. outFile=. 'E:\EMCAS Misc\J Loadflow\lineFLows11.txt' > > > > > > > > rawData=. 1 1000 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 > 0 > > 10 0 11 _1000 > > bus=. 11 2 $ rawData > > > > rawData=. 1 2 4 0.02 2000 1 3 4 0.025 2000 2 > > 3 2 0.08 2000 2 4 3 0.01 2000 2 5 2 0.02 > > 2000 3 8 4 0.04 2000 3 9 2 0.05 2000 4 5 2 > > 0.01 2000 4 6 4 0.02 2000 4 7 3 0.01 2000 5 > > 7 3 0.015 2000 6 7 2 0.01 2000 8 10 4 > > 0.025 2000 8 9 3 0.03 2000 9 10 2 0.04 2000 6 > > 11 3 0.02 2000 7 11 3 0.025 2000 10 11 2 > > 0.04 2000 > > branch=. 18 5 $ rawData > > > > > > netInj=. 1{"1 bus > > > > > > circuits=. 2{"1 branch > > reactance=. 3{"1 branch > > > > > > fromTo=. 0 1{"1 branch NB. select branch from to nodes > > fromTo=. fromTo - 1 NB. subtract 1 for zero based indexing > > fromTo=. ,;/ fromTo > > > > NB. create a nxn matrix for reactance with all zero values > > size =. #bus > > spReactance=. (size,size) $ 0 > > > > > > NB. amend the values in spReactance with values from reactance at > locations > > indicated by fromTo > > spReactance=. reactance fromTo } spReactance > > > > NB. spReactance is upper traingle. transpose and add > > NB. why the hook +|: is not equal to y + (|:y) ? > > tmp=. |: spReactance > > spReactance=. spReactance + tmp > > > > > > NB. create a nxn matrix for no. of circuits > > spCircuits=. (size,size) $ 0 > > spCircuits=. circuits fromTo } spCircuits > > tmp=. |: spCircuits > > spCircuits=. spCircuits + tmp > > > > > > lineReactance=. spReactance%spCircuits > > lineSusceptance=._1 % lineReactance > > > > NB. above two lines can be combined as shown below. But, need > lineReactance > > later > > NB. lineSusceptance=. _1*spCircuits%spReactance > > > > > > NB. replace all __ with 0 > > lineSusceptance=. (lineSusceptance=__)}lineSusceptance,:0 > > > > setDiag =: (>: * i.) @#@] } > > > > NB. line susceptance with diagonal calculations > > lineSusceptance =: (_1 * +/ lineSusceptance) setDiag lineSusceptance > > > > NB. eliminate slack bus and column. Currently last row and last column > > NB. In future should be deleting row and column specified by slack bus > > > > NB. delete last row and column > > lineSusceptance =: _1 }. "1 ] _1 }. lineSusceptance > > > > > > NB. inverse the matrix > > lineSusceptance =: %. lineSusceptance > > > > NB. insert slack row and column with zero values (auto fills in 2nd > > dimension) > > lineSusceptance =: lineSusceptance, ((1, size) $0) > > > > > > NB. define matrix product > > mp=: +/ .* > > > > NB. calculate bus phase angles > > phase =: lineSusceptance mp netInj > > > > > > NB. calculate line flows > > Pij=: -/~@[ % ] > > > > lineFlows=. phase Pij lineReactance > > > > NB. replace all _ and __ with 0 > > lineFlows=. (lineFlows=_)}lineFlows,:0 > > lineFlows=. (lineFlows=__)}lineFlows,:0 > > > > ]lineFlows > > > > NB. write the results to a file > > NB. lineFlows fwrites outFile > > > > > > NB. Check and identify the lines that are over loaded > > NB. TODO: use sparse matrices if possible > > NB. TODO: there are some division by zeros. Can this be avoided? > > NB. TODO: Handle arbitrary bus labels instead of 1 to n > > > > > > On Tue, May 21, 2013 at 1:42 PM, P T <[email protected]> wrote: > > > >> That works perfectly! All the assumptions you made are correct and > thanks > >> for reading my mind. > >> > >> Yes, I am dividing by zero and all those values are not applicable as > >> there are no connections between those nodes in the network. After > cleaning > >> up, I will post the complete program along with sample data. It still > looks > >> like a Java program rather than a J program. > >> > >> Thanks, > >> PT > >> > >> > >> On Tue, May 21, 2013 at 1:11 PM, Raul Miller <[email protected] > >wrote: > >> > >>> On Tue, May 21, 2013 at 12:49 PM, P T <[email protected]> wrote: > >>> > I am almost close to finishing an initial version. But, I don't know > >>> how to > >>> > write a verb P(ij) = (phase(i) - phase(j))/lineReactance(ij) where > >>> > >>> Without looking at your data, this smells like: > >>> > >>> Pij=: -/~@[ % ] > >>> > >>> which you use as: > >>> phase Pij lineReactance > >>> > >>> Here, phase is assumed to be 1 dimensional and lineReactance two > >>> dimensional and square. And phase -:&# lineReactance should always be > >>> true. > >>> > >>> I did not try this on your example data since I do not know what the > >>> result would be. > >>> > >>> I will note that you are dividing by zero, which is worrying. > >>> > >>> -- > >>> Raul > >>> ---------------------------------------------------------------------- > >>> For information about J forums see http://www.jsoftware.com/forums.htm > >>> > >> > >> > > ---------------------------------------------------------------------- > > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
