A Simple Circuit Simulator in Scilab

Version 0.9 [Do not trust anything less than Version 1.0 even if it is made by me.]
Report bugs to dilawar@ee.iitb.ac.in

1. Modified Nodal Analysis

We implemented modified nodal analysis (MNA) in Scilab even though I do not like this method. It is very bad for circuits having more than few thousand nodes. It is implemented mostly for my understanding purpose.

Anyway, please make sure your circuit description file defined in Sec 2.1 does not have current sources in parallel with any of the voltage sources. This is a bug in this implementation and we probably will not remove it in near future unless I get excited about it or you pay me for it [wink, wink].You read about MNA on wikipedia or somewhere else. If the version of this implementation is higher than 0.9, then it is assured that this bug has been removed.

2. Scilab Implementation

2.1 Circuit Description file

  • It is more or less like Spice input file.
  • We do not allow small letters to define any of source e.g. use I instead of i to write down a current source.Use R to write down a resistance instead or r. I do not like using small letters except writing blah but you can modify code to allow them.
  • Use any name after letter I, V and R. You can have same names for different devices. We do not care in thisimplementation. It is recommended that you use IXX, VXX, and RXX format.
  • We only handle V, I and R in this implementation but it can extended to handle all kind of conductanceeasily. May be you’d like to do that exercise.
  • Every line must end with semicolon.
  • Comments starts with #.
  • n0 is always ground node. Other nodes can have any names. It is recommended that you start node name with letter ‘n’.
  • We do not check the topology of the circuit. Give correct net-list.
  • We do not take any responsibility for any damage caused by using this implementation. Nonetheless, mybest wishes.

For example, see the following file 2.1.


# Dilawar - dilawar@ee.iitb.ac.in
# September 18, 2010.
# This is net-list of my circuit. It should be read as spice net-list unless
# otherwise stated explicitly e.g. IXX means current source, VXX is voltage
# source, RXX is resistance and so on.
#
# Only difference here is that we do not allow lower-case letters for circuit
# element names.
# n0 is always the ground node.
IxB n3 n0 2.32;
VyA n1 n0 5;
RyA n1 n2 0.1;
RzB n2 n3 0.1;
RzC n2 n0 0.2;
VyB n3 n0 3;

2.2 Parse this file

Following is the SCILAB function which is used to parse this file. It is adequately commented.

// Dilawar, dilawar@ee.iitb.ac.in
// Sep 18, 2010
// This Scilab function should parse file describing a circuit and should return
// list containing V, I and R.
// Only a circuit which has resistance, voltage
// sources and current sources as elements are dealt with
//. We do not take responsibility for
// any harm caused by using this script. Nonetheless my best wishes.

////////////////////////////////////////////////////////////////////////////////
///
//
// About : Parse a file containing circuit elements and return a list containing
// structure. This structure holds information about the current, voltage, and
// resistance available in the circuit.
//
////////////////////////////////////////////////////////////////////////////////
///
function [voltSource, curSource, valResistance, number] =
parseCircuitFile(myCircuitFilePath)

//Error checking.
[lhs, rhs] = argn(0)
if( rhs == 0) then
    error("Expecting file path.");
end

// Open file.
[fd, err] = mopen(myCircuitFilePath, "r");

// Initialize numbers of entity in file.
noOfcomments = 0;
noOfVoltageSources = 0;
noOfCurrentSources = 0;
noOfResistance = 0;

// Create generic structure to store the data.This is not to be stored. Just an
// overview how it should be used.
stComponent = struct('type', 'V_I_R', 'node1', 'nX', 'node0', 'nY', 'val', 00)

// create lists to store data for each type of element.
curSource= list();
voltSource = list();
valResistance = list();

// Now we parse the file to get the components.
while 0 == meof(fd) then // till end of file is not reached.

    line = mgetl(fd, 1); // get a line.
    if 0 == length(line)
        // empty line, break. Else needed a mechanism to handle this bug.
        break;

    else
        ch = part(line, 1); // get first letter of line.
        if ch == "#" then
            //disp('Comment found.');
            // Do something.
            noOfcomments = noOfcomments + 1;

        elseif ch == "I" then
            //disp('We got a current source.');

            noOfCurrentSources = noOfCurrentSources + 1;

            n = length(line);

            //disp(n);
            // Get the name of the current source.
            countSpace = 0;
            j = 0;
          for i = 1: n,
                //disp('Inside for.', i);
                if part(line, i) == " " then
                    //disp('Found space.');

                    if 0 == countSpace then // first space delimit names.
                        nameCurSource = part(line, 2:i-1);
                        countSpace = countSpace + 1;
                        //disp(nameCurSource);
                        j = i; // remember the location of previous space.

                  elseif    1 == countSpace then // we got a node.
                        upperNode = part(line, j+1:i-1);
                        countSpace = countSpace + 1;
                        j = i;
                        //disp(upperNode);

                    elseif 2 == countSpace then // we got a node
                        lowerNode = part(line, j+1:i-1);
                        countSpace = countSpace + 1;
                        j = i;
                        //disp(lowerNode);
                    end; // inside if-elseif end here.
                end; // outside if ends here.

           // Now there is no space left at the end of line. So we can not read
           // value of the component using this algorithm. We use ';' to read
           // this value.

                if part(line, i) == ";" then // we got a value, uhuu
                    val = part(line, j+1:i-1);
                    countSpace = countSpace + 1;
                    j = i;
                    //disp('Value is:', val);
                end;

            end; // for-loop ends here.

            // Put all the data into matrix.
            stComponent.type = "I"; // this is a current source.
            stComponent.name = nameCurSource;
            stComponent.node1 = upperNode;
            stComponent.node0 = lowerNode;
            stComponent.val = val;

        // This list will contain information in this order. Index [1] stores
        // a matrix in which entries specify that stored is an struct. Index[2]
        // will provide with a matrix [1, 1] which I am not sure what it is.
        // Index[3] is the type of the component. I for current source, V for
        // voltage, R for resistance. Index[4] is name of the component.Index[5]
        // is the upper node of connection of this component in circuit and
        // Index[6] is the lower node. Index[7] is the corresponding value in
        // string which must be converted into double before use.
            curSource = lstcat(curSource, stComponent);

        elseif ch == "R" then
            //disp('We got a resistance.');
            // Do something.
            noOfResistance = noOfResistance + 1;

            n = length(line);

            countSpace = 0;
            j = 0;
          for i = 1: n,
                //disp('Inside for.', i);
                if part(line, i) == " " then
                    //disp('Found space.');

                    if 0 == countSpace then // first space delimit names.
                        nameResistance = part(line, 2:i-1);
                        countSpace = countSpace + 1;
                        //disp(nameResistance);
                        j = i; // remember the location of previous space.

                  elseif    1 == countSpace then // we got a node.
                        upperNode = part(line, j+1:i-1);
                        countSpace = countSpace + 1;
                        j = i;
                        //disp(upperNode);

                    elseif 2 == countSpace then // we got a node
                        lowerNode = part(line, j+1:i-1);
                        countSpace = countSpace + 1;
                        j = i;
                        //disp(lowerNode);
                    end; // inside if-elseif end here.
                end; // outside if ends here.

            // Now there is no space left at the end of line. So we can not read
            // value of the component using this algorithm. We use ';' to read
            // this value.

                if part(line, i) == ";" then // we got a value, uhuu
                    val = part(line, j+1:i-1);
                    countSpace = countSpace + 1;
                    j = i;
                    //disp('Value is:', val);
                end;

            end; // for-loop ends here.

            // Put all the data into matrix.
            stComponent.type = "R"; // this is a resistance.
            stComponent.name = nameResistance;
            stComponent.node1 = upperNode;
            stComponent.node0 = lowerNode;
            stComponent.val = val;

       // This list will contain information in this order. Index 1, stores
       // a matrix in which entries specify that stored is an struct. Index 2
       // will provide with a matrix [1, 1] which I am not sure what it is.
       // Index 3 is the type of the component. I for current source, V for
       // voltage, R for resistance. Index 4  is name of the component.Index 5
       // is the upper node of connection of this component in circuit and
       // Index 6 is the lower node. Index 7 is the corresponding value of the
       // component.
     valResistance = lstcat(valResistance, stComponent);

        elseif ch == "V" then
            //disp('We got a voltage source.');
            // Do something.
            noOfVoltageSources = noOfVoltageSources + 1;

            n = length(line);

            //disp(n);
            // Get the name of the current source.
            countSpace = 0;
            j = 0;
          for i = 1: n,
                //disp('Inside for.', i);
                if part(line, i) == " " then
                    //disp('Found space.');

                    if 0 == countSpace then // first space delimit names.
                        nameVoltSource = part(line, 2:i-1);
                        countSpace = countSpace + 1;
                        //disp(nameVoltSource);
                        j = i; // remember the location of previous space.

                  elseif    1 == countSpace then // we got a node.
                        upperNode = part(line, j+1:i-1);
                        countSpace = countSpace + 1;
                        j = i;
                        //disp(upperNode);

                    elseif 2 == countSpace then // we got a node
                        lowerNode = part(line, j+1:i-1);
                        countSpace = countSpace + 1;
                        j = i;
                        //disp(lowerNode);
                    end; // inside if-elseif end here.
                end; // outside if ends here.

            // Now there is no space left at the end of line. So we can not read
            // value of the component using this algorithm. We use ';' to read
            // this value.

                if part(line, i) == ";" then // we got a value, uhuu
                    val = part(line, j+1:i-1);
                    countSpace = countSpace + 1;
                    j = i;
                    //disp('Value is:', val);
                end;

            end; // for-loop ends here.

            // Put all the data into matrix.
            stComponent.type = "V"; // this is a voltage source.
            stComponent.name = nameVoltSource;
            stComponent.node1 = upperNode;
            stComponent.node0 = lowerNode;
            stComponent.val = val;

          // This list will contain information in this order. Index 1, stores
          // a matrix in which entries specify that stored is an struct. Index 2
          // will provide with a matrix [1, 1] which I am not sure what it is.
          // Index[3] is the type of the component. I for current source, V for
          // voltage, R for resistance. Index 4 is name of the component.Index 5
         // is the upper node of connection of this component in circuit and
         // Index 6 is the lower node. Index 7 is the corresponding value of the
          // component.
            voltSource = lstcat(voltSource, stComponent);

        else
            disp('Mayday. Alien found. Check your circuit file.');
            // Please do something
        end; // if elseif else ends here.
    end; //if else ends here.
end; // while ends here.
mclose(fd);
// return number of elements in circuit.
number = [noOfVoltageSources, noOfCurrentSources, noOfResistance];

endfunction

2.3 Generate coefficient matrix

Well this part should have a better documentation. matM is the matrix which will contain all the coefficients. Matrix matB got all the coefficient matrix. It is better that you work out a simple example for MNA and understand the structure of these matrices.

// Scilab file.
// Funtion
// Dilawar Singh, Sepp 19, 2010.
//
// This script will compute the matrices required to solve the electric circuit
// as described in a text file. To generate the matrix we will use data
// structure returned by a function defined in parseCircuitFile.sce file.
//

// matM : coefficient matrix.
// matB : source matrix.
// listNodes : list having distince nodes at which equations are written.
// nodeV : nodes where voltage sources are attached.
// etc.
function [matM, matB, listNodes, nodeV] = generateMatrix(myFileName)
exec('./parseCircuitFile.sci');
// open file and parse the file.
[vS, iS, rS, numS] = parseCircuitFile(myFileName);
// We got all the sources.

// Build the matrix.
// calculate the no of voltage sources.
nV = numS(1);
// no of current source.
nI = numS(2);
// no of resistors
nR = numS(3);

//////////////////////////////////////////////////
//To do : Extract the topology of the network here.

// First get the node lists.
nodeV = [];
nodeI = [];
nodeR = [];
// node with voltage source.
for j = 1 : nV,
    nodeV = [nodeV, vS((j-1)*7 + 4), vS((j-1)*7 + 5)];
end;
//nodes with current sources.
for j = 1 : nI,
    nodeI = [nodeI, iS((j-1)*7 + 4), iS((j-1)*7 + 5)];
end;
// nodes with resistance
for j = 1 : nR,
    nodeR = [nodeR, rS((j-1)*7 + 4), rS((j-1)*7 + 5)];
end;

// We can now have some sort of circuit topology.

// Extract nodes on which we need to write the nodal equations. This is
// neccessary since a node mostly have repeated entry due to more than one
// element incident to it.
// diffNodes store all such nodes.
// TODO : This is not an efficient implementation.
tempNodes = [nodeV, nodeI, nodeR];
diffNodes = []; // Different node
n = length(length(tempNodes));
i = 1;
allMatch = [];
while i < n then
    j = i +1;
    while j < n then
        if tempNodes(i) ~= tempNodes(j) then
            //disp('Dont Match');
            j = j + 1;
        else
            //disp('Match.');
            allMatch = [allMatch, j];
            tempNodes(j) = []; //Delete it.
            j = j;
            n = length(length(tempNodes));
        end;
    end;
    if length(allMatch) == 0 then
        diffNodes = [diffNodes, tempNodes(i)];
        i = i + 1;
    else
        for p = 1 : length(allMatch),
            tempNodes(j) = [];
        end;
        diffNodes = [diffNodes, tempNodes(i)];
        i = i + 1;
    end;
   n = length(length(tempNodes));
end;

// Discard n0 since by constrain this is our ground node.
// Now we can use any method to solve it. Either MNA or NAL-NBK methods.
// Here we use MNA since it will take less time to implement now.

// Now we will have numNodes + nV variables.
// Add all the resistance incident on one node. They will form the diagonal
// entries for matrix for node equations. Non-diagonal entries can also be
// generated inside this part only. But we give another block of code to make it
// look more straigtforward. Both of these loops should be merged into one for
// less time consumption.
matM = [];
diagM = [];
 // Remove n0 from the diff nodes and use it to build the conducatance matrices.
 // This way it is more natural. We could have ignored this node previously but
 // we keep it, just in case if we need it in future.
 listNodes = [];
 for i = 1 : length(length(diffNodes)),
     if diffNodes(i) ~= 'n0' then
         listNodes = [listNodes, diffNodes(i)];
     end;
 end;

numNodes = length(length(listNodes)); // total nodes excluding ground 'n0'.
for i = 1 : numNodes,
    diagR = 0;
    for p = 1 : (lstsize(rS)/7),
        if listNodes(i) == rS((p-1)*7 + 4) then
            diagR = diagR + 1/evstr(rS((p-1)*7 + 6)); // get the conductance.
        elseif listNodes(i) == rS((p-1)*7 + 5) then
            diagR = diagR + 1/evstr(rS((p-1)*7 + 6));
        end;
    end;
    matM(i, i) = diagR;
end;
// Here we extract non-diagonal entries. Since matM should be symetrical about
// its diagonal. Life is slightly easy.
for i = 1 : numNodes,
    nonDiagR = 0;
    for  j = i+1 : numNodes,
        for p = 1 : lstsize(rS)/7,
            if listNodes(i) == rS((p-1)*7 + 4) then
                if listNodes(j) == rS((p-1)*7 + 5) then
                    nonDiagR = nonDiagR + 1/evstr(rS((p-1)*7 + 6)); //
conductance.
                end;
// it might happen that we have reversed the order of nodes while writing R.
// Though, the way we are extracting the nodes, it should never happen. But
// anyway if someone modifies a part, we should not have any bug here.
            elseif listNodes(i) == rS((p-1)*7 + 5) then
                 if listNodes(j) == rS((p-1)*7 + 4) then
                    // Following is conductance.
                    nonDiagR = nonDiagR + 1/evstr(rS((p-1)*7 +6));
                end;
            end;
        end;
        matM(i, j) = -nonDiagR;
        matM(j, i) = -nonDiagR;
    end;
end;

// Now we need to extend this matM with entries which came from voltage sources.
// Generate a column to be appended to matM for each voltage entry.
// Same time we push these values into matrix B which contains sources. Its
// dimenstion is listNodes+nV by 1.
matB = zeros(numNodes+nV, 1); // mat to hold sources.

for i = 1 : numNodes, // for each node check if there is a voltage source.
    colVoltS = zeros(numNodes,1);
    for j = i+1 : numNodes,
        for p = 1 : lstsize(vS)/7,
            if listNodes(i) == vS((p-1)*7 + 4) then // i is the first node.
                //disp('A voltage source node.');
                if listNodes(j) == vS((p-1)*7 + 5) then // j is the second node.
                    colVoltS(i, 1) = 1;
                    colVoltS(j, 1) = -1;
                    matB(numNodes+i, 1) = evstr(vS((p-1)*7 + 6)); // value V.
                    //disp(colVoltS, 'A');
                else // the second node is n0
                    colVoltS(i, 1) = 1;
                    //disp(colVoltS, 'B');
                    matB(numNodes+i, 1) = evstr(vS((p-1)*7 + 6)); // value V.
                end;
            end;
        end;
    end;
    //disp(colVoltS, 'D');
    // If this vector is non zero. Push it into matM.
    if rank(colVoltS) ~= 0 then
        n = length(matM(1,:));
        for i = 1 : length(colVoltS),
            //disp(length(colVoltS, 'A'));
            matM(i, n + 1) = colVoltS(i);
            matM(n + 1, i) = colVoltS(i);
        end;
    end; // our matrix is ready.
end;

// Final step. We have to push current source into source matrix (vector).
for i = 1 : numNodes, // for each node check if there is a current.
    curS = 0; // to hold current source values.
    for j = i+1 : numNodes,
        for p = 1 : lstsize(iS)/7,
            if listNodes(i) == iS((p-1)*7 + 4) then // i is the first node.
                //disp('A current source node.');
                if listNodes(j) == iS((p-1)*7 + 5) then // j is the second node.
                    curS = curS + evstr(iS((p-1)*7 + 6)); // value i.
                    matB(i, 1) = curS;
                    matB(j, 1) = -curS
                else // the second node is n0
                    curS = curS + evstr(iS((p-1)*7 + 6)); // value i.
                    matB(i, 1) = curS;
                end;
            end;
        end;
    end;
end;

// matM is our coefficient matrix. matB contains sources. Solve it.
msprintf("Your circuit contains %d current sources, %d voltage sources, and ...
%d resistance. PLEASE NOTE THAT at this point we do not check topology of ...
the circuit. Make sure your file describing your circuit is error free.", ...
nI, nV, nR);
endfunction;

2.4 Top most level script and its result

// Scilab file
// Top level script file to solve a circuit using MNA method which I do not like
// very much.

// See respective function files for more details.
// Algorithm implemented is standard MNA method.

// Execute file to load function definition.
exec('./generateMatrix.sci');

// Make sure following file exists.
[matCoeff, matSource, listNodes, nV] = generateMatrix('./myCircuit');
//disp(matCoeff);
//disp(matSource);
//disp(listNodes);
 // solve these equations. matCoeff*x= matSource.
[ansNodes] = linsolve(matCoeff, -matSource);

disp('Node potentials.')
// Display the results.
n = length(length(listNodes));
for i = 1 : n
    mprintf("Node %s votage is : %f\n", listNodes(i), ansNodes(i));
end;
disp('Current through voltage sources.');
for i = n+1 : length(ansNodes),
    mprintf("\nCurrent through voltage source (%s, %s) is : %f", nV(i-n),...
    nV(i-n+1), ansNodes(i));
end;

// This ends our endeavor.
exit();

Results for the file 2.1.

Node potentials.
Node n1 votage is : 5.000000
Node n3 votage is : 3.000000
Node n2 votage is : 3.200000

Current through voltage sources.
Current through voltage source (n1, n0) is : -18.000000
Current through voltage source (n0, n3) is : 4.320000

RESOURCES : mna_scilab_ver0.9

[soucecode language=”matlab”]
Advertisements

About Dilawar

Graduate Student at National Center for Biological Sciences, Bangalore.
This entry was posted in Electrical Network, Technology and Engineering and tagged , , , . Bookmark the permalink.

One Response to A Simple Circuit Simulator in Scilab

  1. Pingback: vincent jackson

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s