C++ Vector4d::bottomRows方法代码示例

  rlc_tank(double r, double l, double c, double vin) : r_(r), l_(l), c_(c), vin_(vin) {
    // state variables are:
    // 0 - input voltage
    // 1 - output voltage
    // 2 - current from input voltage source
    // 3 - current through inductor

    // Create matrices implementing the equation G*X + C*dX/dt = u(t)
    // where u(t) is the independent sources

    Matrix4d C = Matrix4d::Zero(), G = Matrix4d::Zero();

    stamp(G, 0, 1, 1.0/r_);
    stamp(C, 1, c_);
    stamp_i(G, 0, 2);   // connect input current to voltage source
    // for the inductor, two operations:
    stamp(C, 3, l_);    // set derivative coefficient
    stamp_i(G, 1, 3);   // connect inductor current derivative to associated voltage

    // input application matrix
    Matrix<double, 4, 1> B; B << 0, 0, -1, 0;   // match Va up with V0

    // output observation matrix - we want to see input, output, and input current
    Matrix<double, 4, 3> L; L << 1, 0, 0
                               , 0, 1, 0
                               , 0, 0, -1       // extract input current as *out* of the source
                               , 0, 0, 0 ;

    // Rewrite equation as C*dX/dt = -G*X + B*u(t)
    // Further the variables we want to view are Y = L.transpose() * X

    // Now we have a set of differential equations, not all of which are in the required form
    // for ODEINT.  Specifically, some of the derivative terms have a coefficient of zero,
    // so they cannot be meaningfully integrated.  This problem is discussed in detail in:
    // Chen, "A Practical Regularization Technique for Modified Nodal Analysis...",
    // IEEE TCAD, July 2012 and my chosen solution is the simple one used in Su,
    // "Efficient Approximate Balanced Truncation of General Large-Scale RLC Systems via Krylov Methods"
    // Proc. 15th ASP-DAC, 2002

    // Permute C to move the non-zero rows of C to the top while creating a "permutation" matrix
    // we can use to adjust G and B for the state variable reordering.

    // Use Eigen reductions to find zero rows
    auto zero_rows = (C.array() == 0.0).rowwise().all();   // per row "all zeros"

    PermutationMatrix<4, 4, std::size_t> permut;
    std::size_t i, j;
    for (i = 0, j=3; i < j;) {
      // loop invariant: rows > j are all zero; rows < i are not
      while ((i < 4) && !zero_rows(i)) ++i;
      while ((j > 0) && zero_rows(j)) --j;
      if (i < j) {
        // exchange rows i and j via the permutation vector
        permut.applyTranspositionOnTheRight(i, j);
        ++i; --j;
      }
    }

    // 2. Apply permutation to MNA matrices
    Matrix4d Cprime = permut * C * permut;       // permute rows and columns
    Matrix4d Gprime = permut * G * permut;
    Vector4d Bprime = permut * B;    // permute only rows
    Matrix<double, 4, 3> Lprime = permut * L;

    // now the first nonzero_count rows of Cprime, Gprime, and Bprime contain equations acceptable to
    // ODEINT, but the remaining rows do not.  We partition the state space into integrable state (X1)
    // and non-integrable (X2) and do the same with the matrices
    std::size_t zero_count = zero_rows.count();
    std::size_t nonzero_count = 4 - zero_count;

    auto G11 = Gprime.topLeftCorner(nonzero_count, nonzero_count);
    auto G12 = Gprime.topRightCorner(zero_count, zero_count);
    auto G21 = Gprime.bottomLeftCorner(zero_count, zero_count);
    auto G22 = Gprime.bottomRightCorner(zero_count, zero_count);

    auto L1 = Lprime.topRows(nonzero_count);
    auto L2 = Lprime.bottomRows(zero_count);

    auto B1 = Bprime.topRows(nonzero_count);
    auto B2 = Bprime.bottomRows(zero_count);

    // produce reduced equations following Su:
    auto Cred = Cprime.topLeftCorner(nonzero_count, nonzero_count);
    auto G22inv = G22.inverse();   // this is the most expensive thing we will do
    auto Gred = G11 - G12 * G22inv * G21;
    // our L, per PRIMA, is the transpose of the one described in Su
    Lred_ = (L1.transpose() - L2.transpose() * G22inv * G21).transpose();  // simplify?
    auto Bred = B1 - G12 * G22inv * B2;
    // we had no "D" (direct input to output) originally but it can happen, especially
    // since we directly control one state variable that we are mapping to an output
    Dred_ = L2.transpose() * G22inv * B2;

    // Solve to produce a single matrix with the equations for each node
    coeff_ = Cred.ldlt().solve(-1.0 * Gred);

    input_ = Cred.ldlt().solve(Bred);

  }