Merge branch 'develop' into 'master'

add Adams edits in lecture tutorials

See merge request acs/public/teaching/systemtheorie2/lecture-tutorials!4

notebooks/Pysim/Block.py

+import numpy
+sin = numpy.sin
+
+#########################################################################
+#### Abstract class for continuous time block in a block diagram ########
+#########################################################################
+
+class Block(object):
+    def __init__(self):
+        self.ninput = 0    
+        self.noutput = 0    
+        self.nstate = 0
+        self.nparam = 0
+        self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)
+        self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+        self.x = numpy.zeros((self.nstate))
+        self.xo = numpy.zeros((self.nstate))
+        self.u = numpy.zeros((self.ninput))
+        self.y = numpy.zeros((self.noutput))    
+
+    def Init(self,initvalue):
+        self.x = initvalue
+
+    def ChangeParameters(self,value):
+        value = 0
+    
+        
+    def Reset(self):
+        self.x = zeros(self.nstate)
+        
+    def GetInput(self,sc):
+        for i in range(0,self.ninput):
+                self.u[i]=sc.flows[self.inpos[i]]
+    
+    def WriteOutput(self,sc):
+        for i in range(0,self.noutput):
+            sc.flows[self.outpos[i]]=self.y[i]
+            
+    def Step(self,t,dt):
+        for i in range(0,self.nstate):
+                self.x[i] = self.x[i]
+                   
+    def cdxdt(self,x,t):
+        dxdt = numpy.zeros((self.nstate))
+        return dxdt
+        
+    def updatestate(self,t,dt):
+            dx1 = self.cdxdt(self.x,t)
+            xp = numpy.zeros((self.nstate)) 
+            for i in range(0,self.nstate):
+                xp[i] = self.x[i]+dx1[i]*dt      
+            dx2 = self.cdxdt(xp,t+dt)
+            for i in range(0,self.nstate):
+                self.x[i] = self.x[i]+(dx1[i]+dx2[i])/2*dt
+
+            
+            
+#########################################################################
+################### Step function as a signal source ####################
+#########################################################################
+            
+class StepSource(Block):
+    def __init__(self,out,startv,endv,ts):
+            super(StepSource, self).__init__()
+            self.noutput = 1
+            self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+            self.outpos[0] = out
+            self.startv = startv
+            self.endv = endv
+            self.ts = ts
+            self.y = numpy.zeros((self.noutput))   
+            self.y[0] = startv
+                    
+    def Step(self,t,dt):
+        if t>=self.ts:
+            self.y[0] = self.endv
+
+#########################################################################
+####################### Sinusoidal signal source ########################
+#########################################################################
+            
+class SinusoidalSignalSource(Block):
+    def __init__(self,out,startv,Am,om,phi):
+            super(SinusoidalSignalSource, self).__init__()
+            self.noutput = 1
+            self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+            self.outpos[0] = out
+            self.Am = Am
+            self.om = om
+            self.phi = phi
+            self.y = numpy.zeros((self.noutput))   
+            self.y[0] = startv
+                    
+    def Step(self,t,dt):
+            self.y[0] = self.Am*sin(self.om*t+self.phi)    
+            
+
+#########################################################################
+######################### Square signal source ##########################
+#########################################################################
+            
+class SquareSignal(Block):
+    def __init__(self,out,hi,lo,f,duty):
+            super(SquareSignal, self).__init__()
+            self.noutput = 1
+            self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+            self.outpos[0] = out
+            self.hi = hi
+            self.lo = lo
+            self.f = f
+            self.Ts = 1/f
+            self.duty = duty
+            self.y = numpy.zeros((self.noutput))   
+            
+                    
+    def Step(self,t,dt):
+            d = t/self.Ts-floor(t/self.Ts);
+            if d > self.duty:
+                  self.y[0] = self.lo;
+            else:
+                  self.y[0] = self.hi;
+
+
+            
+#########################################################################
+##################### Constant as a signal source #######################
+#########################################################################
+
+class Constant(Block):
+    def __init__(self,out,value):
+            super(Constant, self).__init__()
+            self.noutput = 1
+            self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+            self.outpos[0] = out
+            self.y = numpy.zeros((self.noutput))   
+            self.y[0] = value
+
+
+#########################################################################
+########################## Division operation  ##########################
+#########################################################################
+            
+class Division(Block):
+    def __init__(self,in1,in2,out):
+            super(Division, self).__init__()
+            self.noutput = 1
+            self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+            self.outpos[0] = out
+            self.ninput = 2
+            self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)   
+            self.inpos[0] = in1
+            self.inpos[1] = in2
+            self.y[0] = 0
+                    
+    def Step(self,t,dt):
+        self.y[0] = self.u[0]/self.u[1]
+
+        
+        
+#########################################################################
+########################## Product operation  ###########################
+#########################################################################
+
+class Product(Block):
+    def __init__(self,in1,in2,out):
+            super(Product, self).__init__()
+            self.noutput = 1
+            self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+            self.outpos[0] = out
+            self.ninput = 2
+            self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)   
+            self.inpos[0] = in1
+            self.inpos[1] = in2
+            self.y = numpy.zeros((self.noutput))   
+            self.y[0] = 0
+                    
+    def Step(self,t,dt):
+        self.y[0] = self.u[0]*self.u[1]
+
+        
+#########################################################################
+############################# Sum operation  ############################
+#########################################################################
+
+class Sum(Block):
+    def __init__(self,in1,in2,out,sign1,sign2):
+            super(Sum, self).__init__()
+            self.noutput = 1
+            self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+            self.outpos[0] = out
+            self.ninput = 2
+            self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)   
+            self.inpos[0] = in1
+            self.inpos[1] = in2
+            self.sign1 = sign1
+            self.sign2 = sign2
+            self.u = numpy.zeros((self.ninput)) 
+            self.y = numpy.zeros((self.noutput))   
+            self.y[0] = 0
+                    
+    def Step(self,t,dt):
+        self.y[0] = self.u[0]*self.sign1+self.u[1]*self.sign2
+
+#########################################################################
+############################## Saturation  ##############################
+#########################################################################
+
+class Saturation(Block):
+    def __init__(self,inp,out,minout,maxout):
+            super(Saturation, self).__init__()
+            self.noutput = 1
+            self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+            self.outpos[0] = out
+            self.ninput = 1
+            self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)   
+            self.inpos[0] = inp
+            self.minout = minout
+            self.maxout = maxout
+            self.u = numpy.zeros((self.ninput)) 
+            self.y = numpy.zeros((self.noutput))   
+
+                    
+    def Step(self,t,dt):
+        self.y[0] = self.u[0]
+        if self.y[0] > self.maxout:
+            self.y[0] = self.maxout
+        if self.y[0] < self.minout:
+            self.y[0] = self.minout
+     
+    def ChangeParameters(self,value):
+        self.minout = value[0]
+        self.maxout = value[1]
+        
+
+#########################################################################
+############################# Gain block  ###############################
+#########################################################################
+
+class Gain(Block):
+    def __init__(self,inp,out,k):
+            super(Gain, self).__init__()
+            self.ninput = len(inp)
+            self.noutput = len(out)
+            self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+            self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)   
+            self.u = numpy.zeros((self.ninput))
+            self.y = numpy.zeros((self.noutput)) 
+            for i in range(0,self.ninput):
+                self.inpos[i] = inp[i]
+            for i in range(0,self.noutput):
+                self.outpos[i] = out[i]   
+            self.k = numpy.zeros((self.noutput,self.ninput))
+            self.k = k
+    
+    def Step(self,t,dt):
+        self.y = self.k.dot(self.u.transpose())
+        
+    def ChangeParameters(self,value):
+        self.k = value
+
+
+#########################################################################
+########################## PI Controller  ###############################
+#########################################################################
+
+class PI(Block):
+    def __init__(self,inp,out,Kp,Ki,ini):                
+        super(PI, self).__init__()
+        self.ninput = 1
+        self.noutput = 1
+        self.nstate = 1
+        self.A = 0
+        self.B = 1
+        self.C = Ki
+        self.D = Kp 
+        self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)
+        self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+        self.x = numpy.zeros((self.nstate))
+        self.xo = numpy.zeros((self.nstate))
+        self.u = numpy.zeros((self.ninput))
+        self.y = numpy.zeros((self.noutput)) 
+        self.inpos[0] = inp
+        self.outpos[0] = out
+        self.x[0] = ini
+        self.xo[0] = ini
+              
+    def Step(self,t,dt):
+        self.y[0] = self.C*self.x[0]+self.D*self.u[0]
+        self.updatestate(t,dt)                
+        
+    def Reset(self):
+        self.x[0] = xo[0]
+                                   
+    def cdxdt(self,x,t):    
+        dxdt = numpy.zeros((self.nstate))
+        dxdt[0] = self.B*u[0]
+        return dxdt
+    
+    def ChangeParameters(self,value):
+        self.C = np.array([value[1]])
+        self.D = np.array([value[0]])
+
+
+#########################################################################
+########################## PID Controller  ##############################
+#########################################################################
+
+class PID(Block):
+    def __init__(self,inp,out,Kp,Ki,Kd,ini):                
+        super(PID, self).__init__()
+        self.ninput = 1
+        self.noutput = 1
+        self.nstate = 1
+        self.B = 1
+        self.C = Ki
+        self.D = Kp 
+        self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)
+        self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+        self.x = numpy.zeros((self.nstate))
+        self.xo = numpy.zeros((self.nstate))
+        self.u = numpy.zeros((self.ninput))
+        self.y = numpy.zeros((self.noutput)) 
+        self.inpos[0] = inp
+        self.outpos[0] = out
+        self.x[0] = ini
+        self.xo[0] = ini
+        self.olderr = 0
+        self.Kd = Kd
+              
+    def Step(self,t,dt):
+        self.y[0] = self.C*self.x[0]+self.D*self.u[0]
+        self.y[0] = self.y[0]+self.Kd*(self.u[0]-self.olderr)/dt
+        self.olderr = self.u[0]
+        self.updatestate(t,dt)                
+        
+    def Reset(self):
+        self.x[0] = xo[0]
+        self.olderr = 0;
+                                   
+    def cdxdt(self,x,t):
+        dxdt = numpy.zeros((self.nstate))
+        dxdt[0] = self.B*self.u[0]
+        return dxdt
+    
+    def ChangeParameters(self,value):
+        self.C = value[1]
+        self.D = value[0]
+        self.Kd = value[2]
+        
+        
+
+#########################################################################
+############################# Integrator ################################
+#########################################################################
+
+class Integrator(Block):
+    def __init__(self,inp,out,ini):                
+        super(Integrator, self).__init__()
+        self.ninput = 1
+        self.noutput = 1
+        self.nstate = 1
+        self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)
+        self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+        self.x = numpy.zeros((self.nstate))
+        self.xo = numpy.zeros((self.nstate))
+        self.u = numpy.zeros((self.ninput))
+        self.y = numpy.zeros((self.noutput)) 
+        self.inpos[0] = inp
+        self.outpos[0] = out
+        self.x[0] = ini
+        self.xo[0] = ini
+                      
+    def Step(self,t,dt):
+        self.y[0] = self.x[0]
+        self.updatestate(t,dt)                
+        
+    def Reset(self):
+        self.x[0] = xo[0]
+                                   
+    def cdxdt(self,x,t):
+        dxdt = numpy.zeros((self.nstate))
+        dxdt[0] = self.u[0]
+        return dxdt
+    
+
+#########################################################################
+################################ Not ####################################
+#########################################################################
+
+class Not(Block):
+    def __init__(self,inp,out):                
+        super(Not, self).__init__()
+        self.ninput = 1
+        self.noutput = 1
+        self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)
+        self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+        self.u = numpy.zeros((self.ninput))
+        self.y = numpy.zeros((self.noutput)) 
+        self.inpos[0] = inp
+        self.outpos[0] = out
+                              
+    def Step(self,t,dt):
+        if self.u[0] == 0:
+            self.y[0] = 1
+        else:    
+            self.y[0] = 0
+        
+
+#########################################################################
+########################## State Space  #################################
+#########################################################################
+
+class StateSpace(Block):
+    def __init__(self,inp,out,A,B,C,D,xo):                
+        super(StateSpace, self).__init__()
+        self.ninput = len(inp)
+        self.noutput = len(out)
+        self.nstate = len(A)
+        self.A = numpy.zeros((self.nstate,self.nstate))
+        self.B = numpy.zeros((self.nstate,self.ninput))
+        self.C = numpy.zeros((self.noutput,self.nstate))
+        self.D = numpy.zeros((self.noutput,self.ninput)) 
+        self.A = A
+        self.B = B
+        self.C = C
+        self.D = D
+        self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)
+        self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+        self.x = numpy.zeros((self.nstate))
+        self.xo = numpy.zeros((self.nstate))
+        self.u = numpy.zeros((self.ninput))
+        self.y = numpy.zeros((self.noutput)) 
+        for i in range(0,self.ninput):
+            self.inpos[i] = inp[i]
+        for i in range(0,self.noutput):
+            self.outpos[i] = out[i]
+        for i in range(0,self.nstate):
+            self.x[i] = xo[i]
+            self.xo[i] = xo[i]
+       
+        
+    def Step(self,t,dt):
+        for i in range(0,self.noutput):
+            self.y[i] = 0
+            for j in range(0,self.nstate):
+                self.y[i] = self.y[i] + self.C[j]*self.x[j]
+            for j in range(0,self.ninput):
+                self.y[i] = self.y[i] +self.D[j]*self.u[j]
+        self.updatestate(t,dt)                
+        
+    def Reset(self):
+        for i in range(0,self.nstate):
+            self.x[i] = xo[i]
+                                   
+    def cdxdt(self,x,t):     
+        dxdt = numpy.zeros((self.nstate)) 
+        for i in range (0,self.nstate):
+            for j in range (0,self.ninput):
+                dxdt[i] = self.B[i,j]*self.u[j]
+            for j in range (0,self.nstate):
+                dxdt[i] = dxdt[i] + self.A[i,j]*x[j]
+        return dxdt
+        
+
+#########################################################################
+####################### Transfer Function  ##############################
+#########################################################################
+
+class TransferFunction(Block):
+    def __init__(self,inp,out,num,den):                
+        super(TransferFunction, self).__init__()
+        self.ninput = 1
+        self.noutput = 1
+        self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)
+        self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+        self.inpos[0] = inp
+        self.outpos[0] = out
+        self.nstate = len(den)-1
+        denord = len(den)
+        numord = len(num)
+        self.A = numpy.zeros((self.nstate,self.nstate))
+        self.B = numpy.zeros((self.nstate))
+        self.C = numpy.zeros((self.nstate))
+        self.D = 0 
+        scale = den[0];
+        for i in range(0,denord):
+            den[i]= den[i]/scale
+        for i in range(0,numord):
+            num[i]= num[i]/scale
+        if numord==denord:
+            self.D = num[0]
+            for i in range(0,numord-1):
+                self.C[i]=num[numord-i-1]-num[0]*den[numord-i-1]
+        else:
+            for i in range(0,numord):
+                self.C[i]=num[numord-i-1]
+        for i in range(0,self.nstate):
+            self.A[self.nstate-1,i] = -den[denord-i-1]
+        for i in range(0,self.nstate-1):
+            self.A[i,i+1]=1
+        self.B[self.nstate-1] = 1
+        self.x = numpy.zeros((self.nstate))
+        self.xo = numpy.zeros((self.nstate))
+        self.u = numpy.zeros((self.ninput))
+        self.y = numpy.zeros((self.noutput)) 
+        for i in range(0,self.nstate):
+            self.x[i] = 0
+            self.xo[i] = 0     
+        
+    def Step(self,t,dt):
+        self.y[0] = self.D*self.u[0]
+        for i in range (0,self.nstate):
+            self.y[0] = self.y[0]+self.C[i]*self.x[i]
+        self.updatestate(t,dt)                
+        
+    def Reset(self):
+        for i in range(0,self.nstate):
+            self.x[i] = xo[i]
+                                   
+    def cdxdt(self,x,t):     
+        dxdt = numpy.zeros((self.nstate)) 
+        for i in range (0,self.nstate):
+            dxdt[i] = self.B[i]*self.u[0]
+            for j in range (0,self.nstate):
+                dxdt[i] = dxdt[i] + self.A[i,j]*x[j]
+        return dxdt
+        
+        

notebooks/Pysim/Circuit.py

+# -*- coding: utf-8 -*-
+import numpy
+
+class LinProb(object):
+    def __init__(self, nodeNum):        
+        self.G = numpy.zeros((nodeNum,nodeNum), dtype=numpy.complex)
+        self.b = numpy.zeros((nodeNum,1), dtype=numpy.complex)
+
+    def getSize(self):
+        return self.G.shape[0]
+        
+    def resetbVector(self):        
+        self.b = numpy.zeros((self.b.shape[0],1), dtype=numpy.complex)
+
+    def solve(self):
+        result = numpy.dot(numpy.linalg.inv(self.G), self.b)
+        return result
+        
+    def show(self):
+        print(self.G)
+        print(self.b)
+
+class Circuit(object):
+    def __init__(self, nnode):
+        self.n = nnode      
+        self.NetList = list()
+        
+    def addComponent(self,h):
+        print("adding component" + str(h))
+        self.nelements = len(self.NetList)+1
+        self.NetList.append(h)
+        
+    def initialize(self,dt):
+        self.Pn = LinProb(self.n)
+        self.Pn.show()
+        self.nelements = len(self.NetList)        
+        print("number of components " + str(self.nelements))
+        for i in range(0,self.nelements):
+            print("initializing component " + str(i+1))
+            self.NetList[i].initialize(self.Pn,dt)
+            
+        return self.Pn
+    
+    def initializeDp(self,dt,om):
+        self.Pn = LinProb(self.n)
+        self.nelements = len(self.NetList)
+        for i in range(0,self.nelements):
+            self.NetList[i].initialize(self.Pn,dt,om)
+            
+    def solveLinTrans(self,tini,tfin,dt):
+        print("show problem")
+        self.Pn.show()        
+        t = tini
+        k = 0
+        
+        while t < tfin:
+            #print("solving for t = " + str(t))
+            self.Pn.resetbVector()            
+            
+            for i in range(0,self.nelements):
+                self.Pn = self.NetList[i].step(self.Pn,dt,t)
+                
+            if t == 0:
+                self.vout = self.Pn.solve()
+            else:
+                self.vout = numpy.append(self.vout, self.Pn.solve(), 1)            
+            
+            for i in range(0,self.nelements):
+                self.NetList[i].postStep(self.vout[:,k],dt)
+                
+            t = t+dt
+            k = k+1
+            
+        return self.vout
+    
+        
+        
+        
+        
\ No newline at end of file

notebooks/Pysim/DTBlock.py

+import numpy
+from Block import Block
+
+#########################################################################
+#### Abstract class for discrete time block in a block diagram ########
+#########################################################################
+
+class DTBlock(Block):
+    def __init__(self):
+        self.ninput = 0    
+        self.noutput = 0    
+        self.nstate = 0
+        self.nparam = 0
+        self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)
+        self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+        self.x = numpy.zeros((self.nstate))
+        self.xo = numpy.zeros((self.nstate))
+        self.u = numpy.zeros((self.ninput))
+        self.y = numpy.zeros((self.noutput))    
+        self.Ts = 0.001
+        self.lastt = -10
+        
+    def Step(self,t,dt):
+        if t-self.lastt >= self.Ts:
+            self.updatediscrete()
+            self.lastt = t
+                     
+    def updatediscrete(self):
+        self.y[0] = 0
+        
+    def Reset(self):
+        self.lastt = -10;
+
+#########################################################################
+################################ FIR ####################################
+#########################################################################
+
+class FIR(DTBlock):
+    def __init__(self,inp,out,a,ts):
+            super(FIR, self).__init__()
+            self.ninput = 1
+            self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)
+            self.inpos[0] = inp            
+            self.noutput = 1
+            self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+            self.outpos[0] = out
+            self.Ts = ts
+            self.pt = len(a)
+            self.a = numpy.zeros((self.pt))
+            self.a = a
+            self.buf = numpy.zeros((self.pt))  
+            self.u = numpy.zeros((self.ninput)) 
+            self.y = numpy.zeros((self.noutput))   
+            self.y[0] = 0
+                    
+    def updatediscrete(self):
+            for i in range (self.pt-1,0,-1):
+                self.buf[i] = self.buf[i-1]
+            self.buf[0] = self.u[0]
+            self.y[0] = 0
+            for i in range (0,self.pt):
+                self.y[0] = self.y[0]+self.a[i]*self.buf[i]      
+           
+    def Reset(self):
+        self.buf = zeros(self.pt,1)
+        self.lastt = -10
+            
+
+
+#########################################################################
+#########################Discrete Integral###############################
+#########################################################################
+
+class DTIntegral(DTBlock):
+    def __init__(self,inp,out,ts,xo):
+            super(DTIntegral, self).__init__()
+            self.ninput = 1
+            self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)
+            self.inpos[0] = inp            
+            self.noutput = 1
+            self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+            self.outpos[0] = out
+            self.Ts = ts
+            self.xo = xo
+            self.x = xo
+            self.u = numpy.zeros((self.ninput)) 
+            self.y = numpy.zeros((self.noutput))   
+            self.y[0] = xo
+                    
+    def updatediscrete(self):
+           self.y[0] = self.x
+           self.x = self.x + self.u[0]*self.Ts
+           
+    def Reset(self):
+        self.x = self.xo  
+
+
+
+#########################################################################
+#########################   ZOH         #################################
+#########################################################################
+
+class ZOH(DTBlock):
+    def __init__(self,inp,out,ts):
+            super(ZOH, self).__init__()
+            self.ninput = 1
+            self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)
+            self.inpos[0] = inp            
+            self.noutput = 1
+            self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+            self.outpos[0] = out
+            self.Ts = ts
+            self.u = numpy.zeros((self.ninput)) 
+            self.y = numpy.zeros((self.noutput))   
+            self.y[0] = 0
+                    
+    def updatediscrete(self):
+           self.y[0] = self.u[0]
+           
+
+
+        
+#########################################################################
+#######################   Delay         #################################
+#########################################################################
+
+class DTDelay(DTBlock):
+    def __init__(self,inp,out,init,ts):
+        super(DTDelay, self).__init__()
+        self.ninput = 1
+        self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)
+        self.inpos[0] = inp            
+        self.noutput = 1
+        self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+        self.outpos[0] = out
+        self.Ts = ts
+        self.u = numpy.zeros((self.ninput)) 
+        self.y = numpy.zeros((self.noutput))   
+        self.y[0] = 0
+        self.oldv = init
+        self.init = init
+                    
+    def updatediscrete(self):
+        self.y[0] = self.oldv
+        self.oldv = self.u[0]
+
+    def ChangeParameters(self,value):
+        self.Ts = value
+            
+    def Reset(self):
+        self.oldv = self.init
+        
+
+        
+        
+        
+#########################################################################
+################################ IIR ####################################
+#########################################################################
+
+class IIR(DTBlock):
+    def __init__(self,inp,out,a,b,ts):
+            super(IIR, self).__init__()
+            self.ninput = 1
+            self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)
+            self.inpos[0] = inp            
+            self.noutput = 1
+            self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+            self.outpos[0] = out
+            self.Ts = ts
+            self.pn = len(a)
+            self.a = numpy.zeros((self.pn))
+            self.a = a
+            self.pd = len(b)
+            self.b = numpy.zeros((self.pd))
+            self.b = b
+            self.bufn = numpy.zeros((self.pn))
+            self.bufd = numpy.zeros((self.pd))
+            self.u = numpy.zeros((self.ninput)) 
+            self.y = numpy.zeros((self.noutput))   
+            self.y[0] = 0
+                    
+    def updatediscrete(self):
+            for i in range (self.pn-1,0,-1):
+                self.bufn[i] = self.bufn[i-1]
+            self.bufn[0] = self.u[0]
+            self.y[0] = 0
+            for i in range (0,self.pn):
+                self.y[0] = self.y[0]+self.a[i]*self.bufn[i]               
+            for i in range (self.pd-1,0,-1):
+                self.bufd[i] = self.bufd[i-1]
+            self.bufd[0] = self.y[0]
+            for i in range (0,self.pd):
+                self.y[0] = self.y[0]-self.b[i]*self.bufd[i]      
+          
+            
+    def Reset(self):
+        self.bufn = zeros(self.pn,1)
+        self.bufd = zeros(self.pd,1)
+        self.lastt = -10
+            
+
+
+
+#########################################################################
+############### Discrete Time Transfer Function  ########################
+#########################################################################
+
+class DTTransferFunction(DTBlock):
+    def __init__(self,inp,out,num,den,Ts):                
+        super(DTTransferFunction, self).__init__()
+        self.Ts = Ts
+        self.ninput = 1
+        self.noutput = 1
+        self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)
+        self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+        self.inpos[0] = inp
+        self.outpos[0] = out
+        self.nstate = len(den)-1
+        denord = len(den)
+        numord = len(num)
+        self.A = numpy.zeros((self.nstate,self.nstate))
+        self.B = numpy.zeros((self.nstate))
+        self.C = numpy.zeros((self.nstate))
+        self.D = 0 
+        scale = den[0];
+        for i in range(0,denord):
+            den[i]= den[i]/scale
+        for i in range(0,numord):
+            num[i]= num[i]/scale
+        if numord==denord:
+            self.D = num[0]
+            for i in range(0,numord-1):
+                self.C[i]=num[numord-i-1]-num[0]*den[numord-i-1]
+        else:
+            for i in range(0,numord):
+                self.C[i]=num[numord-i-1]
+        for i in range(0,self.nstate):
+            self.A[self.nstate-1,i] = -den[denord-i-1]
+        for i in range(0,self.nstate-1):
+            self.A[i,i+1]=1
+        self.B[self.nstate-1] = 1
+        self.x = numpy.zeros((self.nstate))
+        self.u = numpy.zeros((self.ninput))
+        self.y = numpy.zeros((self.noutput)) 
+        for i in range(0,self.nstate):
+            self.x[i] = 0  
+        
+    def updatediscrete(self):
+        self.y[0] = self.D*self.u[0]
+        for i in range (0,self.nstate):
+            self.y[0] = self.y[0]+self.C[i]*self.x[i]
+        xn = numpy.zeros((self.nstate)) 
+        for i in range (0,self.nstate):
+            xn[i] = self.B[i]*self.u[0]
+            for j in range (0,self.nstate):
+                xn[i] = xn[i] + self.A[i,j]*self.x[j]
+        for i in range (0,self.nstate):
+            self.x[i] = xn[i]
+            
+        
+    def Reset(self):
+        for i in range(0,self.nstate):
+            self.x[i] = 0
+        
+
+
+#########################################################################
+########################### DT State Space  #############################
+#########################################################################
+
+class DTStateSpace(DTBlock):
+    def __init__(self,inp,out,A,B,C,D,xo,Ts):                
+        super(DTStateSpace, self).__init__()
+        self.Ts = Ts
+        self.ninput = len(inp)
+        self.noutput = len(out)
+        self.nstate = len(A)
+        self.inpos = numpy.zeros((self.ninput),dtype=numpy.int)
+        self.outpos = numpy.zeros((self.noutput),dtype=numpy.int)
+        self.A = numpy.zeros((self.nstate,self.nstate))
+        self.B = numpy.zeros((self.nstate,self.ninput))
+        self.C = numpy.zeros((self.nstate,self.noutput))
+        self.D = numpy.zeros((self.noutput,self.ninput))
+        self.x = numpy.zeros((self.nstate))
+        self.u = numpy.zeros((self.ninput))
+        self.y = numpy.zeros((self.noutput)) 
+        for i in range(0,self.ninput):
+            self.inpos[i] = inp[i]
+        for i in range(0,self.noutput):
+            self.outpos[i] = out[i]
+        for i in range(0,self.nstate):
+            self.x[i] = xo[i]
+            self.xo[i] = xo[i]
+        
+    def updatediscrete(self):
+        for i in range(0,self.noutput):
+            self.y[i] = 0
+            for j in range(0,self.nstate):
+                self.y[i] = self.y[i] + self.C[j]*self.x[j]
+            for j in range(0,self.ninput):
+                self.y[i] = self.y[i] +self.D[j]*self.u[j]
+        xn = numpy.zeros((self.nstate)) 
+        for i in range (0,self.nstate):
+            for j in range (0,self.ninput):
+                xn[i] = xn[i]+self.B[i,j]*self.u[j]
+            for j in range (0,self.nstate):
+                xn[i] = xn[i] + self.A[i,j]*self.x[j]
+        for i in range (0,self.nstate):
+            self.x[i] = xn[i]
+            
+        
+    def Reset(self):
+        for i in range(0,self.nstate):
+            self.x[i] = 0
+        
+        

notebooks/Pysim/RCComponents.py

+# -*- coding: utf-8 -*-
+
+import numpy
+
+class RCComponent(object):
+    def __init__(self, p, n):
+        self.Pos = p
+        self.Neg = n
+        self.linear = 1
+        self.hyin = 0
+        self.hyout = 0
+        
+    def applyMatrixStamp(self, P,dt):
+        return P
+        
+    def initialize(self, P,dt): 
+        return P     
+        
+    def step(self,P,dt,t):
+        return P 
+
+    def postStep(self,vsol,dt):
+        return vsol
+
+class Resistance(RCComponent):           
+    def __init__(self,p,n,par):      
+        super(Resistance, self).__init__(p,n)
+        self.r = par        
+        if par > 0:
+            self.Gc = 1/par
+        else:
+            print("Error: Resistance value not correct")
+            
+    def initialize(self,P,dt):
+        print("initializing resistance")
+        return Resistance.applyMatrixStamp(self,P,dt)
+            
+    def applyMatrixStamp(self,P,dt):  
+        if self.Pos > -1:
+            P.G[self.Pos,self.Pos] = P.G[self.Pos,self.Pos] + self.Gc
+        if self.Neg > -1:
+            P.G[self.Neg,self.Neg] = P.G[self.Neg,self.Neg] + self.Gc
+        if (self.Pos > -1) and (self.Neg > -1):
+            P.G[self.Pos,self.Neg] = P.G[self.Pos,self.Neg] - self.Gc
+            P.G[self.Neg,self.Pos] = P.G[self.Neg,self.Pos] - self.Gc  
+        P.show()              
+        return P         
+        
+            
+class Capacitor(RCComponent):
+    def __init__(self,p,n,par):
+        super(Capacitor, self).__init__(p,n)
+        self.C = par
+        self.Gl = 0
+        
+    def initialize(self,P,dt):
+        print("initializing capacitor")
+        self.Bc = 0
+        self.Vc = 0
+        self.Ic = 0
+        return Capacitor.applyMatrixStamp(self,P,dt)
+        
+    def applyMatrixStamp(self,P,dt): 
+        self.Gc = 2*self.C/dt      
+        if self.Pos > -1:
+            P.G[self.Pos,self.Pos] = P.G[self.Pos,self.Pos] + self.Gc
+        if self.Neg > -1:
+            P.G[self.Neg,self.Neg] = P.G[self.Neg,self.Neg] + self.Gc
+        if (self.Pos > -1) and (self.Neg > -1):
+            P.G[self.Pos,self.Neg] = P.G[self.Pos,self.Neg] - self.Gc
+            P.G[self.Neg,self.Pos] = P.G[self.Neg,self.Pos] - self.Gc      
+            
+        P.show()   
+        return P      
+        
+    def step(self,P,dt,t):
+        self.Bc = self.Gc*self.Vc+self.Ic
+        if self.Pos > -1:
+            P.b[self.Pos] = P.b[self.Pos] + self.Bc
+        if self.Neg > -1:
+            P.b[self.Neg] = P.b[self.Neg] - self.Bc
+        return P
+        
+    def postStep(self,vsol,dt):
+        if self.Pos > -1:
+            if self.Neg > -1:
+                self.Vc = vsol[self.Pos] - vsol[self.Neg]
+            else:
+                self.Vc = vsol[self.Pos]              
+        else:
+            self.Vc = -vsol[self.Neg];    
+                    
+        self.Ic = self.Gc*self.Vc-self.Bc;
+
+        
+class Inductor(RCComponent):    
+    def __init__(self,p,n,par):
+        super(Inductor, self).__init__(p,n)
+        self.L = par
+        
+    def initialize(self,P,dt):        
+        print("initializing inductor")
+        self.Bl = 0
+        self.Vl = 0
+        self.Il = 0
+        return Inductor.applyMatrixStamp(self,P,dt)
+        
+    def applyMatrixStamp(self,P,dt): 
+        self.Gl = dt/(2*self.L)      
+        if self.Pos > -1:
+            P.G[self.Pos,self.Pos] = P.G[self.Pos,self.Pos] + self.Gl
+        if self.Neg > -1:
+            P.G[self.Neg,self.Neg] = P.G[self.Neg,self.Neg] + self.Gl
+        if (self.Pos > -1) and (self.Neg > -1):
+            P.G[self.Pos,self.Neg] = P.G[self.Pos,self.Neg] - self.Gl
+            P.G[self.Neg,self.Pos] = P.G[self.Neg,self.Pos] - self.Gl   
+        P.show()             
+        return P       
+        
+    def step(self,P,dt,t):
+        self.Bl = self.Gl*self.Vl+self.Il
+        if self.Pos > -1:
+            P.b[self.Pos] = P.b[self.Pos] - self.Bl
+        if self.Neg > -1:
+            P.b[self.Neg] = P.b[self.Neg] + self.Bl
+        return P
+                
+    def postStep(self,vsol,dt):
+        if self.Pos > -1:
+            if self.Neg > -1:
+                self.Vl = vsol[self.Pos] - vsol[self.Neg]
+            else:
+                self.Vl = vsol[self.Pos]              
+        else:
+            self.Vl = -vsol[self.Neg];    
+                    
+        self.Il = self.Gl*self.Vl+self.Bl;
+        
+class IdealACVoltageSource(RCComponent):    
+    def __init__(self,p,n,V,omega,phase):
+        super(IdealACVoltageSource, self).__init__(p,n)
+        self.Vs = V
+        self.om = omega
+        self.ph = phase
+        
+    def initialize(self,P,dt):
+        print("initializing ideal ac voltage source")
+        return IdealACVoltageSource.applyMatrixStamp(self,P,dt)
+        
+    def applyMatrixStamp(self,P,dt):         
+        self.ntot = P.getSize()
+        self.possource = self.ntot
+        print("position source at " + str(self.possource))
+        P.G = numpy.append(P.G, numpy.zeros((self.ntot,1), dtype=numpy.complex), 1)
+        P.G = numpy.append(P.G, numpy.zeros((1,self.ntot+1), dtype=numpy.complex), 0)
+        P.b = numpy.append(P.b, numpy.zeros((1,1), dtype=numpy.complex), 0)
+                
+        if self.Pos > -1:
+            P.G[self.possource,self.Pos] = 1
+            P.G[self.Pos,self.possource] = 1
+        if self.Neg > -1:
+            P.G[self.possource,self.Neg] = -1
+            P.G[self.Neg,self.possource] = -1  
+            
+        P.show()
+        return P       
+            
+    def step(self,P,dt,t):
+        P.b[self.possource,0] = self.Vs * numpy.sin(self.om*t+self.ph)
+        return P
+        
\ No newline at end of file

notebooks/Pysim/Schema.py

+# -*- coding: utf-8 -*-
+import numpy as np
+
+
+class Schema(object):
+    def __init__(self, ti,tf,dt,numflows):
+        self.tini = ti    
+        self.tfin = tf    
+        self.nmodels = 0
+        self.nflows = numflows
+        self.dt = dt
+        self.flows = np.zeros(numflows)
+        self.BlockList = list()
+        
+    def AddComponent(self,h):
+        self.nmodels = len(self.BlockList)+1
+        self.BlockList.append(h)
+
+    def AddListComponents(self,l):
+        p = len(l)
+        for i in range(0,p):
+            self.BlockList.append(l[i])
+            self.nmodels = self.nmodels+1
+
+    def Init(self):
+        self.t = self.tini
+        for i in range(0,self.nflows):
+            self.flows[i] = 0
+
+    def Reset(self):
+        for i in range(0,self.nmodels):
+            self.BlockList[i].Reset()
+
+    def Run(self,listout):
+        self.Init()
+        qmax = listout.size            
+        npoints = (self.tfin-self.tini)/self.dt
+        npt = int(npoints)
+        out = np.zeros((qmax+1,npt-1))
+        p=0
+        while self.Step() and p<npt-1:
+            out[0,p] = self.GetTime()
+            for q in range(0,qmax):
+                out[q+1, p] = self.GetFlow(listout[q])
+            p=p+1
+        return out
+            
+    def ChangeTimeStep(self,dt):
+        self.dt = dt
+
+    def ChangeParameters(self,pos,val):
+        self.ModelList[pos].ChangeParameters(val)
+        
+    def Step(self):
+        for  i in range(0,self.nmodels): 
+            self.BlockList[i].GetInput(self)
+            self.BlockList[i].Step(self.t,self.dt)
+            self.BlockList[i].WriteOutput(self)
+        self.t = self.t+self.dt
+        if self.t < self.tfin:
+                flag = 1
+        else:
+                flag = 0
+        return flag    
+    
+    def GetFlow(self,n):
+        val = self.flows[n]
+        return val     
+           
+    def GetTime(self):
+        val = self.t
+        return val
+            
\ No newline at end of file

notebooks/V07_steuerbarkeit_und_beobachtbarkeit.ipynb

@@ -7,6 +7,32 @@
     "# <span style='color:OrangeRed'>V7 - Steuerbarkeit und Beobachtbarkeit</span>"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "pi = np.pi\n",
+    "sin = np.sin\n",
+    "cos = np.cos\n",
+    "sqrt = np.sqrt\n",
+    "exp = np.exp\n",
+    "atan2 = np.arctan2\n",
+    "log10 = np.log10\n",
+    "\n",
+    "def print_det(M):\n",
+    "    print(\"Determinante: \"+str(np.around(np.linalg.det(M))))\n",
+    "\n",
+    "def print_rank(M):\n",
+    "    print(\"Rang: \"+str(np.linalg.matrix_rank(S_C)))"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -37,11 +63,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 30,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
-    "A = [-1 -1 -2; 0 -1 1; 1 0 -1];"
+    "A = np.array([[-1,-1,-2],[0,-1,1],[1,0,-1]]);"
    ]
   },
   {
@@ -58,11 +84,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 31,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [],
    "source": [
-    "b = [2; 0; 1];"
+    "b = np.array([[2],[0],[1]]);"
    ]
   },
   {
@@ -77,11 +103,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 32,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [],
    "source": [
-    "c = [1 1 0];"
+    "c = np.array([1,1,0]);"
    ]
   },
   {
@@ -125,58 +151,60 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 33,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "S_C =\n",
-      "\n",
-      "   2  -4   1\n",
-      "   0   1   0\n",
-      "   1   1  -5\n",
-      "\n"
+      "S_C:\n",
+      "[[ 2 -4  1]\n",
+      " [ 0  1  0]\n",
+      " [ 1  1 -5]]\n"
      ]
     }
    ],
    "source": [
-    "S_C = [b A*b A*A*b]"
+    "column1 = b\n",
+    "column2 = np.matmul(A,b)\n",
+    "column3 = np.matmul(np.matmul(A,A),b)\n",
+    "S_C = np.c_[column1,column2,column3]\n",
+    "print(\"S_C:\\n\"+str(S_C))"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 34,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "ans = -11\n"
+      "Determinante: -11.0\n"
      ]
     }
    ],
    "source": [
-    "det(S_C)"
+    "print_det(S_C)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 35,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "ans =  3\n"
+      "Rang: 3\n"
      ]
     }
    ],
    "source": [
-    "rank(S_C)"
+    "print_rank(S_C)"
    ]
   },
   {
@@ -212,58 +240,57 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 36,
+   "execution_count": 8,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "S_O =\n",
-      "\n",
-      "   1   1   0\n",
-      "  -1  -2  -1\n",
-      "   0   3   1\n",
-      "\n"
+      "S_C:\n",
+      "[[ 1 -1  0]\n",
+      " [ 1 -2  3]\n",
+      " [ 0 -1  1]]\n"
      ]
     }
    ],
    "source": [
-    "S_O = [c; c*A; c*A*A]"
+    "S_O = np.c_[c, np.matmul(c,A), np.matmul(np.matmul(c,A),A)]\n",
+    "print(\"S_C:\\n\"+str(S_O))"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 37,
+   "execution_count": 9,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "ans =  2\n"
+      "Determinante: 2.0\n"
      ]
     }
    ],
    "source": [
-    "det(S_O)"
+    "print_det(S_O)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 38,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "ans =  3\n"
+      "Rang: 3\n"
      ]
     }
    ],
    "source": [
-    "rank(S_O)"
+    "print_rank(S_O)"
    ]
   },
   {
@@ -274,33 +301,32 @@
     " <code>det(S_O)</code> $\\ne$ 0 und <code>rank(S_O)</code> = $n$ = $3$\n",
     "<br>Damit ist das gegebene System  <span style='color:red'>vollständig beobachtbar</span>! "
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Octave",
-   "language": "octave",
-   "name": "octave"
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
   },
   "language_info": {
-   "file_extension": ".m",
-   "help_links": [
-    {
-     "text": "GNU Octave",
-     "url": "https://www.gnu.org/software/octave/support.html"
-    },
-    {
-     "text": "Octave Kernel",
-     "url": "https://github.com/Calysto/octave_kernel"
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
    },
-    {
-     "text": "MetaKernel Magics",
-     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
-    }
-   ],
-   "mimetype": "text/x-octave",
-   "name": "octave",
-   "version": "4.2.2"
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.6"
   }
  },
  "nbformat": 4,

requirements.txt

 scilab_kernel
 numpy
 scipy
+control
 matplotlib
 PyGnuplot