# Branched System 3-Reservoir Example

Consider the branched system shown in below.  In this example the friction factor is assumed constant for simplicity, but in practice would vary during the solution computations.



The hydraulics question is what is the discharge in each pipe and what is the total head at the junction (notice we don't know the junction elevation in this example --- if the elevation were specified, we could also find the pressure head).

## Solution Approach

First populate the four equations with the appropriate numerical values.

\begin{equation}
70 = h_D + (0.015)\frac{5000}{0.6}\frac{V_{AD}|V_{AD}|}{2(9.8)}
\end{equation}

\begin{equation}
100 = h_D + (0.015)\frac{3000}{0.8}\frac{V_{BD}|V_{BD}|}{2(9.8)}
\end{equation}

\begin{equation}
-80 = - h_D  + (0.015)\frac{4000}{1.2}\frac{V_{DC}|V_{DC}|}{2(9.8)}
\end{equation}

\begin{equation}
\frac{\pi(0.6)^2}{4} V_{AD}+\frac{\pi(0.8)^2}{4} V_{BD} = \frac{\pi(1.2)^2}{4} V_{DC}
\end{equation}

Next, compute all the constants, and organize the 4 equations into a system of simultaneous equations

\begin{equation}
\begin{matrix}
h_D & 6.377 V_{AD}|V_{AD}| & 0 & 0  & = 70\\
h_D & 0 & 2.869 V_{BD}|V_{BD}| & 0  & = 100\\
h_D & 0 & 0 & - 2.551 V_{DC}|V_{DC}|  & = 80\\
0 & 0.2827 V_{AD} & 0.5026 V_{BD}& -1.1309 V_{DC}  & = 0\\
\end{matrix}
\end{equation}

Because the system is non-linear we need to employ some appropriate method, herein we will use a Quasi-Newton method with numerical approximations to derivatives.

First a more conventional vector-matrix structure

\begin{gather}
\begin{pmatrix}
1 & 6.377 |V_{AD}| & 0 & 0 \\      
1 & 0 & 2.869 |V_{BD}| & 0 \\    
1 & 0 & 0 & - 2.551 |V_{DC}| \\    
0 & 0.2827  & 0.5026 & -1.1309 \\ 
\end{pmatrix}
\bullet
\begin{pmatrix}
 h_D \\
V_{AD} \\
V_{BD} \\
V_{DC} \\
\end{pmatrix}
-
\begin{pmatrix}
70 \\
100 \\
 80 \\
0 \\
\end{pmatrix}
=
\begin{pmatrix}
0 \\
0 \\
0 \\
0 \\
\end{pmatrix}
\end{gather}


## Apply the multiple-variable Extension of Newton’s Method

The solution can use the methods directly from ENGR-1330 and similar courses, here we just implement the methods.  A very brief explaination is in the readings, but multi-variable Newton's method is usually taught in calculus.

Here we will define some support modules `linearsolver()` and `vector_matrix_lib` which are just python source code containing prototype functions which are embedded into the notebook, but ordinarily would save them in separate file and use `import`.


```python
# SolveLinearSystem.py
# Code to read A and b
# Then solve Ax = b for x by Gaussian elimination with back substitution
##########
def linearsolver(A,b):
    n = len(A)
# M = A  #this is object to object equivalence
# copy A into M element by element
    M=[[0.0 for jcol in range(n)]for irow in range(n)]
    for irow in range(n):
        for jcol in range(n):
            M[irow][jcol]=A[irow][jcol]
    i = 0
    for x in M:
     x.append(b[i])
     i += 1
    for k in range(n):
     for i in range(k,n):
       if abs(M[i][k]) > abs(M[k][k]):
          M[k], M[i] = M[i],M[k]
       else:
          pass
     for j in range(k+1,n):
         q = float(M[j][k]) / M[k][k]
         for m in range(k, n+1):
            M[j][m] -=  q * M[k][m]
    x = [0 for i in range(n)]
    x[n-1] =float(M[n-1][n])/M[n-1][n-1]
    for i in range (n-1,-1,-1):
      z = 0
      for j in range(i+1,n):
          z = z  + float(M[i][j])*x[j]
      x[i] = float(M[i][n] - z)/M[i][i]
#    print (x)
    return(x)
```


```python
# vector_matrix_lib.py
# useful linear algebra tools
import math   # This will import math module

def writeM(M,ir,jc,label):
    print ("------",label,"------")
    for i in range(0,ir,1):
        print (M[i][0:jc])
    print ("-----------------------------")
    #return

def writeV(V,ir,label):
    print ("------",label,"------")
    for i in range(0,ir,1):
        print (V[i])
    print ("-----------------------------")
    #return

def mmmult(amatrix,bmatrix,rowNumA,colNumA,rowNumB,colNumB):
    AB =[[0.0 for j in range(colNumB)] for i in range(rowNumA)]
    for i in range(0,rowNumA):
        for j in range(0,colNumB):
            for k in range(0,colNumA):
                AB[i][j]=AB[i][j]+amatrix[i][k]*bmatrix[k][j]
    return(AB)

def mvmult(amatrix,xvector,rowNumA,colNumA):
    bvector=[0.0 for i in range(rowNumA)]
    for i in range(0,rowNumA):
        for j in range(0,1):
            for k in range(0,colNumA):
                bvector[i]=bvector[i]+amatrix[i][k]*xvector[k]
    return(bvector)

def vvadd(avector,bvector,length):
    aplusb=[0.0 for i in range(length)]
    for i in range(length):
        aplusb[i] = avector[i] + bvector[i]
    return(aplusb)

def vvsub(avector,bvector,length):
    aminusb=[0.0 for i in range(length)]
    for i in range(length):
        aminusb[i] = avector[i] - bvector[i]
    return(aminusb)
             
def vdotv(avector,bvector,length):
    adotb=0.0
    for i in range(length):
        adotb=adotb+avector[i]*bvector[i]
    return(adotb)
```



Next we define the four equations, as well a a vector valued function `func` and its Jacobian `jacob` as below.  


```python
#################################################################
# Newton Solver Example -- Numerical Derivatives                #
#################################################################
import math   # This will import math module from python distribution
#from LinearSolverPivot import linearsolver # This will import our solver module
#from vector_matrix_lib import writeM,writeV,vdotv,vvsub # This will import our vector functions

def eq1(x1,x2,x3,x4):
    eq1 = 1*x1 + 6.377*x2 +0*x3 + 0*x4 - 70
    return(eq1)
def eq2(x1,x2,x3,x4):
    eq2 = 1*x1 + 0*x2 +2.869*x3 + 0*x4 - 100
    return(eq2)
def eq3(x1,x2,x3,x4):
    eq3 = 1*x1 + 0*x2 +0*x3 + -2.551*x4 -80
    return(eq3)
def eq4(x1,x2,x3,x4):
    eq4 = 0*x1 + 0.2827*x2 +0.5026*x3 + -1.1309*x4 - 0
    return(eq4)
##############################################################
def func(x1,x2,x3,x4):
    func  = [0.0 for i in range(4)] # null list
    # build the function
    func[0] = eq1(x1,x2,x3,x4)
    func[1] = eq2(x1,x2,x3,x4)
    func[2] = eq3(x1,x2,x3,x4)
    func[3] = eq4(x1,x2,x3,x4)
    return(func)

def jacob(x1,x2,x3,x4):
    jacob = [[0.0 for j in range(4)] for i in range(4)] # constructed list  
    #build the  jacobian
    jacob[0][0]=1
    jacob[0][1]=6.377*x2
    jacob[0][2]=0
    jacob[0][3]=0
    jacob[1][0]=1
    jacob[1][1]=0
    jacob[1][2]=2.869*x3
    jacob[1][3]=0
    jacob[2][0]=1
    jacob[2][1]=0
    jacob[2][2]=0
    jacob[2][3]=-2.551*x4
    jacob[3][0]=0
    jacob[3][1]=0.2827
    jacob[3][2]=0.5026
    jacob[3][3]=-1.1309
    return(jacob)
```

Next we create vectors to store values, and supply initial guesses to the system, and echo the inputs.


```python
deltax  = [0.0 for i in range(4)] # null list
xguess  = [0.0 for i in range(4)] # null list
myfunc  = [0.0 for i in range(4)] # null list
myjacob = [[0.0 for j in range(4)] for i in range(4)] # constructed list  
# supply initial guess
xguess[0] = float(input("Value for x1 : "))
xguess[1] = float(input("Value for x2 : "))
xguess[2] = float(input("Value for x3 : "))
xguess[3] = float(input("Value for x4 : "))
# build the initial function
myfunc = func(xguess[0],xguess[1],xguess[2],xguess[3])
#build the initial jacobian
myjacob=jacob(xguess[0],xguess[1],xguess[2],xguess[3])
#write initial results
writeV(xguess,4,"Initial X vector ")
writeV(myfunc,4,"Initial FUNC vector ")
writeM(myjacob,4,4,"Initial Jacobian ")
# solver parameters
tolerancef = 1.0e-9
tolerancex = 1.0e-9
# Newton-Raphson
```

    Value for x1 :  1
    Value for x2 :  1
    Value for x3 :  1
    Value for x4 :  1


    ------ Initial X vector  ------
    1.0
    1.0
    1.0
    1.0
    -----------------------------
    ------ Initial FUNC vector  ------
    -62.623
    -96.131
    -81.551
    -0.3455999999999999
    -----------------------------
    ------ Initial Jacobian  ------
    [1, 6.377, 0, 0]
    [1, 0, 2.869, 0]
    [1, 0, 0, -2.551]
    [0, 0.2827, 0.5026, -1.1309]
    -----------------------------


Now we apply the algorithm a few times, here the count is set to 10. So eneter the loop, test for stopping, then update.


```python
for iteration in range(24):
    myfunc = func(xguess[0],xguess[1],xguess[2],xguess[3])
    testf = vdotv(myfunc,myfunc,4)
    if testf <= tolerancef :
        print("f(x) close to zero\n test value : ", testf)
        break
    myjacob=jacob(xguess[0],xguess[1],xguess[2],xguess[3])
    deltax=linearsolver(myjacob,myfunc)
    testx = vdotv(deltax,deltax,4)
    if testx <= tolerancex :
        print("solution change small\n test value : ", testx)
        break
    xguess=vvsub(xguess,deltax,4)
##    print("iteration : ",iteration)
##    writeV(xguess,2,"Current X vector ")
##    writeV(myfunc,2,"Current FUNC vector ")
print("Exiting Iteration : ",iteration)
writeV(xguess,4,"Exiting X vector ")
writeV(myfunc,4,"Exiting FUNC vector using Finite-Differences")
```

    f(x) close to zero
     test value :  2.0273725264180003e-28
    Exiting Iteration :  1
    ------ Exiting X vector  ------
    84.61708956275575
    -2.292157685864162
    5.361767318663032
    1.8099135879089594
    -----------------------------
    ------ Exiting FUNC vector using Finite-Differences ------
    -1.4210854715202004e-14
    0.0
    0.0
    -8.881784197001252e-16
    -----------------------------


Now tidy up the output, convert back to discharges (from velocities) and report results.


```python
print("Head at Junction : %.3f" % xguess[0], " feet ")
print("Discharge in Pipe AD: %.3f" % (0.2827 * xguess[1]), " cubic feet per second")
print("Discharge in Pipe BD: %.3f" % (0.5026 * xguess[2]), " cubic feet per second ")
print("Discharge in Pipe DC: %.3f" % (1.1309 * xguess[3]), " cubic feet per second ")
```

    Head at Junction : 84.617  feet 
    Discharge in Pipe AD: -0.648  cubic feet per second
    Discharge in Pipe BD: 2.695  cubic feet per second 
    Discharge in Pipe DC: 2.047  cubic feet per second 



```python

```

00ceed7cce53826061f5e9fba971b11d49a63a6e
ipynb
Jupyter Notebook
9-MyJupyterNotebooks/15-NonLinearSystems/BranchedPipelineExample.ipynb
dustykat/engr-1330-psuedo-course
3e7e31a32a1896fcb1fd82b573daa5248e465a36
CC0-1.0
9-MyJupyterNotebooks/15-NonLinearSystems/BranchedPipelineExample.ipynb
dustykat/engr-1330-psuedo-course
3e7e31a32a1896fcb1fd82b573daa5248e465a36
CC0-1.0
9-MyJupyterNotebooks/15-NonLinearSystems/BranchedPipelineExample.ipynb
dustykat/engr-1330-psuedo-course
3e7e31a32a1896fcb1fd82b573daa5248e465a36
CC0-1.0
Qwen/Qwen-72B
1. YES
2. YES