Allokering utav minne - Acceptabelt inom inbyggda system?

[DanielM]

Jag håller på med MATLAB till C-kodgenerering. Det fungerar riktigt fint och ger en datorgenererad kod. Men koden är riktigt komplex och innehåller mycket free, malloc och calloc. Är detta acceptabelt inom inbyggda system, eller bör jag undvika det så långt som det går? Standarden är C89/C90.

Notera att C-kodgenerering genererar väldigt mycket kod och skulle man skriva koden själv så skulle det bli en mycket mindre kod.
Detta är bara main filen.

Kod: Markera allt

#include "rt_nonfinite.h"
#include "mpc.h"
#include "main.h"
#include "mpc_terminate.h"
#include "mpc_emxAPI.h"
#include "mpc_initialize.h"

/* Function Declarations */
static emxArray_real_T *argInit_Unboundedx1_real_T(void);
static double argInit_real_T(void);
static emxArray_real_T *c_argInit_UnboundedxUnbounded_r(void);
static void main_mpc(void);

/* Function Definitions */
static emxArray_real_T *argInit_Unboundedx1_real_T(void)
{
  emxArray_real_T *result;
  static int iv1[1] = { 2 };

  int idx0;

  /* Set the size of the array.
     Change this size to the value that the application requires. */
  result = emxCreateND_real_T(1, iv1);

  /* Loop over the array to initialize each element. */
  for (idx0 = 0; idx0 < result->size[0U]; idx0++) {
    /* Set the value of the array element.
       Change this value to the value that the application requires. */
    result->data[idx0] = argInit_real_T();
  }

  return result;
}

static double argInit_real_T(void)
{
  return 0.0;
}

static emxArray_real_T *c_argInit_UnboundedxUnbounded_r(void)
{
  emxArray_real_T *result;
  static int iv0[2] = { 2, 2 };

  int idx0;
  int idx1;

  /* Set the size of the array.
     Change this size to the value that the application requires. */
  result = emxCreateND_real_T(2, iv0);

  /* Loop over the array to initialize each element. */
  for (idx0 = 0; idx0 < result->size[0U]; idx0++) {
    for (idx1 = 0; idx1 < result->size[1U]; idx1++) {
      /* Set the value of the array element.
         Change this value to the value that the application requires. */
      result->data[idx0 + result->size[0] * idx1] = argInit_real_T();
    }
  }

  return result;
}

static void main_mpc(void)
{
  emxArray_real_T *U;
  emxArray_real_T *A;
  emxArray_real_T *B;
  emxArray_real_T *C;
  emxArray_real_T *x;
  emxInitArray_real_T(&U, 1);

  /* Initialize function 'mpc' input arguments. */
  /* Initialize function input argument 'A'. */
  A = c_argInit_UnboundedxUnbounded_r();

  /* Initialize function input argument 'B'. */
  B = c_argInit_UnboundedxUnbounded_r();

  /* Initialize function input argument 'C'. */
  C = c_argInit_UnboundedxUnbounded_r();

  /* Initialize function input argument 'x'. */
  x = argInit_Unboundedx1_real_T();

  /* Call the entry-point 'mpc'. */
  mpc(A, B, C, x, argInit_real_T(), argInit_real_T(), argInit_real_T(), U);
  emxDestroyArray_real_T(U);
  emxDestroyArray_real_T(x);
  emxDestroyArray_real_T(C);
  emxDestroyArray_real_T(B);
  emxDestroyArray_real_T(A);
}

int main(int argc, const char * const argv[])
{
  (void)argc;
  (void)argv;

  /* Initialize the application.
     You do not need to do this more than one time. */
  mpc_initialize();

  /* Invoke the entry-point functions.
     You can call entry-point functions multiple times. */
  main_mpc();

  /* Terminate the application.
     You do not need to do this more than one time. */
  mpc_terminate();
  return 0;
}

/* End of code generation (main.c) */

Huvudfilen ser ut så här.

Kod: Markera allt

/* Include files */
#include "rt_nonfinite.h"
#include "mpc.h"
#include "mpc_emxutil.h"
#include "mpower.h"

/* Function Declarations */
static void gammaMat(const emxArray_real_T *A, const emxArray_real_T *B, const
                     emxArray_real_T *C, double N, emxArray_real_T *GAMMA);

/* Function Definitions */
static void gammaMat(const emxArray_real_T *A, const emxArray_real_T *B, const
                     emxArray_real_T *C, double N, emxArray_real_T *GAMMA)
{
  int i1;
  int loop_ub;
  int i;
  emxArray_real_T *r1;
  emxArray_real_T *r2;
  emxArray_real_T *y;
  emxArray_real_T *CAB;
  emxArray_real_T *varargin_1;
  int m;
  int inner;
  int coffset;
  int i2;
  int boffset;
  int b_i;
  int c_i;
  int k;
  int aoffset;
  double d0;
  double b_N;
  boolean_T empty_non_axis_sizes;

  /*  Create the lower triangular toeplitz matrix */
  i1 = GAMMA->size[0] * GAMMA->size[1];
  GAMMA->size[0] = (int)((double)C->size[0] * N);
  GAMMA->size[1] = (int)((double)B->size[1] * N);
  emxEnsureCapacity_real_T(GAMMA, i1);
  loop_ub = (int)((double)C->size[0] * N) * (int)((double)B->size[1] * N);
  for (i1 = 0; i1 < loop_ub; i1++) {
    GAMMA->data[i1] = 0.0;
  }

  i = 0;
  emxInit_real_T(&r1, 2);
  emxInit_real_T(&r2, 2);
  emxInit_real_T(&y, 2);
  emxInit_real_T(&CAB, 2);
  emxInit_real_T(&varargin_1, 2);
  while (i <= (int)N - 1) {
    if ((C->size[1] == 1) || (A->size[0] == 1)) {
      i1 = y->size[0] * y->size[1];
      y->size[0] = C->size[0];
      y->size[1] = A->size[1];
      emxEnsureCapacity_real_T(y, i1);
      loop_ub = C->size[0];
      for (i1 = 0; i1 < loop_ub; i1++) {
        coffset = A->size[1];
        for (i2 = 0; i2 < coffset; i2++) {
          y->data[i1 + y->size[0] * i2] = 0.0;
          boffset = C->size[1];
          for (c_i = 0; c_i < boffset; c_i++) {
            y->data[i1 + y->size[0] * i2] += C->data[i1 + C->size[0] * c_i] *
              A->data[c_i + A->size[0] * i2];
          }
        }
      }
    } else {
      m = C->size[0];
      inner = C->size[1];
      i1 = y->size[0] * y->size[1];
      y->size[0] = C->size[0];
      y->size[1] = A->size[1];
      emxEnsureCapacity_real_T(y, i1);
      for (loop_ub = 0; loop_ub < A->size[1]; loop_ub++) {
        coffset = loop_ub * m - 1;
        boffset = loop_ub * inner - 1;
        for (b_i = 1; b_i <= m; b_i++) {
          y->data[coffset + b_i] = 0.0;
        }

        for (k = 1; k <= inner; k++) {
          if (A->data[boffset + k] != 0.0) {
            aoffset = (k - 1) * m;
            for (b_i = 1; b_i <= m; b_i++) {
              y->data[coffset + b_i] += A->data[boffset + k] * C->data[(aoffset
                + b_i) - 1];
            }
          }
        }
      }
    }

    if ((y->size[1] == 1) || (B->size[0] == 1)) {
      i1 = r1->size[0] * r1->size[1];
      r1->size[0] = y->size[0];
      r1->size[1] = B->size[1];
      emxEnsureCapacity_real_T(r1, i1);
      loop_ub = y->size[0];
      for (i1 = 0; i1 < loop_ub; i1++) {
        coffset = B->size[1];
        for (i2 = 0; i2 < coffset; i2++) {
          r1->data[i1 + r1->size[0] * i2] = 0.0;
          boffset = y->size[1];
          for (c_i = 0; c_i < boffset; c_i++) {
            r1->data[i1 + r1->size[0] * i2] += y->data[i1 + y->size[0] * c_i] *
              B->data[c_i + B->size[0] * i2];
          }
        }
      }
    } else {
      m = y->size[0];
      inner = y->size[1];
      i1 = r1->size[0] * r1->size[1];
      r1->size[0] = y->size[0];
      r1->size[1] = B->size[1];
      emxEnsureCapacity_real_T(r1, i1);
      for (loop_ub = 0; loop_ub < B->size[1]; loop_ub++) {
        coffset = loop_ub * m - 1;
        boffset = loop_ub * inner - 1;
        for (b_i = 1; b_i <= m; b_i++) {
          r1->data[coffset + b_i] = 0.0;
        }

        for (k = 1; k <= inner; k++) {
          if (B->data[boffset + k] != 0.0) {
            aoffset = (k - 1) * m;
            for (b_i = 1; b_i <= m; b_i++) {
              r1->data[coffset + b_i] += B->data[boffset + k] * y->data[(aoffset
                + b_i) - 1];
            }
          }
        }
      }
    }

    if ((C->size[1] == 1) || (A->size[0] == 1)) {
      i1 = y->size[0] * y->size[1];
      y->size[0] = C->size[0];
      y->size[1] = A->size[1];
      emxEnsureCapacity_real_T(y, i1);
      loop_ub = C->size[0];
      for (i1 = 0; i1 < loop_ub; i1++) {
        coffset = A->size[1];
        for (i2 = 0; i2 < coffset; i2++) {
          y->data[i1 + y->size[0] * i2] = 0.0;
          boffset = C->size[1];
          for (c_i = 0; c_i < boffset; c_i++) {
            y->data[i1 + y->size[0] * i2] += C->data[i1 + C->size[0] * c_i] *
              A->data[c_i + A->size[0] * i2];
          }
        }
      }
    } else {
      m = C->size[0];
      inner = C->size[1];
      i1 = y->size[0] * y->size[1];
      y->size[0] = C->size[0];
      y->size[1] = A->size[1];
      emxEnsureCapacity_real_T(y, i1);
      for (loop_ub = 0; loop_ub < A->size[1]; loop_ub++) {
        coffset = loop_ub * m - 1;
        boffset = loop_ub * inner - 1;
        for (b_i = 1; b_i <= m; b_i++) {
          y->data[coffset + b_i] = 0.0;
        }

        for (k = 1; k <= inner; k++) {
          if (A->data[boffset + k] != 0.0) {
            aoffset = (k - 1) * m;
            for (b_i = 1; b_i <= m; b_i++) {
              y->data[coffset + b_i] += A->data[boffset + k] * C->data[(aoffset
                + b_i) - 1];
            }
          }
        }
      }
    }

    if ((y->size[1] == 1) || (B->size[0] == 1)) {
      i1 = r2->size[0] * r2->size[1];
      r2->size[0] = y->size[0];
      r2->size[1] = B->size[1];
      emxEnsureCapacity_real_T(r2, i1);
      loop_ub = y->size[0];
      for (i1 = 0; i1 < loop_ub; i1++) {
        coffset = B->size[1];
        for (i2 = 0; i2 < coffset; i2++) {
          r2->data[i1 + r2->size[0] * i2] = 0.0;
          boffset = y->size[1];
          for (c_i = 0; c_i < boffset; c_i++) {
            r2->data[i1 + r2->size[0] * i2] += y->data[i1 + y->size[0] * c_i] *
              B->data[c_i + B->size[0] * i2];
          }
        }
      }
    } else {
      m = y->size[0];
      inner = y->size[1];
      i1 = r2->size[0] * r2->size[1];
      r2->size[0] = y->size[0];
      r2->size[1] = B->size[1];
      emxEnsureCapacity_real_T(r2, i1);
      for (loop_ub = 0; loop_ub < B->size[1]; loop_ub++) {
        coffset = loop_ub * m - 1;
        boffset = loop_ub * inner - 1;
        for (b_i = 1; b_i <= m; b_i++) {
          r2->data[coffset + b_i] = 0.0;
        }

        for (k = 1; k <= inner; k++) {
          if (B->data[boffset + k] != 0.0) {
            aoffset = (k - 1) * m;
            for (b_i = 1; b_i <= m; b_i++) {
              r2->data[coffset + b_i] += B->data[boffset + k] * y->data[(aoffset
                + b_i) - 1];
            }
          }
        }
      }
    }

    d0 = (double)C->size[0] * (1.0 + (double)i);
    if (1.0 + (double)i > d0) {
      i1 = 1;
    } else {
      i1 = i + 1;
    }

    b_N = (N - (1.0 + (double)i)) + 1.0;

    /*  Create the column for the GAMMA matrix */
    i2 = CAB->size[0] * CAB->size[1];
    CAB->size[0] = (int)((double)C->size[0] * b_N);
    CAB->size[1] = B->size[1];
    emxEnsureCapacity_real_T(CAB, i2);
    loop_ub = (int)((double)C->size[0] * b_N) * B->size[1];
    for (i2 = 0; i2 < loop_ub; i2++) {
      CAB->data[i2] = 0.0;
    }

    for (b_i = 0; b_i < (int)((b_N - 1.0) + 1.0); b_i++) {
      d0 = (double)C->size[0] * ((double)b_i + 1.0);
      if ((double)b_i + 1.0 > d0) {
        i2 = 1;
      } else {
        i2 = b_i + 1;
      }

      mpower(A, b_i, varargin_1);
      if ((C->size[1] == 1) || (varargin_1->size[0] == 1)) {
        c_i = y->size[0] * y->size[1];
        y->size[0] = C->size[0];
        y->size[1] = varargin_1->size[1];
        emxEnsureCapacity_real_T(y, c_i);
        loop_ub = C->size[0];
        for (c_i = 0; c_i < loop_ub; c_i++) {
          coffset = varargin_1->size[1];
          for (k = 0; k < coffset; k++) {
            y->data[c_i + y->size[0] * k] = 0.0;
            boffset = C->size[1];
            for (aoffset = 0; aoffset < boffset; aoffset++) {
              y->data[c_i + y->size[0] * k] += C->data[c_i + C->size[0] *
                aoffset] * varargin_1->data[aoffset + varargin_1->size[0] * k];
            }
          }
        }
      } else {
        m = C->size[0];
        inner = C->size[1];
        c_i = y->size[0] * y->size[1];
        y->size[0] = C->size[0];
        y->size[1] = varargin_1->size[1];
        emxEnsureCapacity_real_T(y, c_i);
        for (loop_ub = 0; loop_ub < varargin_1->size[1]; loop_ub++) {
          coffset = loop_ub * m - 1;
          boffset = loop_ub * inner - 1;
          for (c_i = 1; c_i <= m; c_i++) {
            y->data[coffset + c_i] = 0.0;
          }

          for (k = 1; k <= inner; k++) {
            if (varargin_1->data[boffset + k] != 0.0) {
              aoffset = (k - 1) * m;
              for (c_i = 1; c_i <= m; c_i++) {
                y->data[coffset + c_i] += varargin_1->data[boffset + k] *
                  C->data[(aoffset + c_i) - 1];
              }
            }
          }
        }
      }

      if ((y->size[1] == 1) || (B->size[0] == 1)) {
        c_i = varargin_1->size[0] * varargin_1->size[1];
        varargin_1->size[0] = y->size[0];
        varargin_1->size[1] = B->size[1];
        emxEnsureCapacity_real_T(varargin_1, c_i);
        loop_ub = y->size[0];
        for (c_i = 0; c_i < loop_ub; c_i++) {
          coffset = B->size[1];
          for (k = 0; k < coffset; k++) {
            varargin_1->data[c_i + varargin_1->size[0] * k] = 0.0;
            boffset = y->size[1];
            for (aoffset = 0; aoffset < boffset; aoffset++) {
              varargin_1->data[c_i + varargin_1->size[0] * k] += y->data[c_i +
                y->size[0] * aoffset] * B->data[aoffset + B->size[0] * k];
            }
          }
        }
      } else {
        m = y->size[0];
        inner = y->size[1];
        c_i = varargin_1->size[0] * varargin_1->size[1];
        varargin_1->size[0] = y->size[0];
        varargin_1->size[1] = B->size[1];
        emxEnsureCapacity_real_T(varargin_1, c_i);
        for (loop_ub = 0; loop_ub < B->size[1]; loop_ub++) {
          coffset = loop_ub * m - 1;
          boffset = loop_ub * inner - 1;
          for (c_i = 1; c_i <= m; c_i++) {
            varargin_1->data[coffset + c_i] = 0.0;
          }

          for (k = 1; k <= inner; k++) {
            if (B->data[boffset + k] != 0.0) {
              aoffset = (k - 1) * m;
              for (c_i = 1; c_i <= m; c_i++) {
                varargin_1->data[coffset + c_i] += B->data[boffset + k] *
                  y->data[(aoffset + c_i) - 1];
              }
            }
          }
        }
      }

      loop_ub = varargin_1->size[1];
      for (c_i = 0; c_i < loop_ub; c_i++) {
        coffset = varargin_1->size[0];
        for (k = 0; k < coffset; k++) {
          CAB->data[((i2 + k) + CAB->size[0] * c_i) - 1] = varargin_1->data[k +
            varargin_1->size[0] * c_i];
        }
      }
    }

    i2 = varargin_1->size[0] * varargin_1->size[1];
    varargin_1->size[0] = (int)(((1.0 + (double)i) - 1.0) * (double)r1->size[0]);
    varargin_1->size[1] = r2->size[1];
    emxEnsureCapacity_real_T(varargin_1, i2);
    loop_ub = (int)(((1.0 + (double)i) - 1.0) * (double)r1->size[0]) * r2->size
      [1];
    for (i2 = 0; i2 < loop_ub; i2++) {
      varargin_1->data[i2] = 0.0;
    }

    coffset = (int)(((1.0 + (double)i) - 1.0) * (double)r1->size[0]);
    boffset = r2->size[1];
    if (!((coffset == 0) || (boffset == 0))) {
      boffset = r2->size[1];
    } else if (!((CAB->size[0] == 0) || (CAB->size[1] == 0))) {
      boffset = CAB->size[1];
    } else {
      boffset = r2->size[1];
      if (boffset > 0) {
        boffset = r2->size[1];
      } else {
        boffset = 0;
      }

      if (CAB->size[1] > boffset) {
        boffset = CAB->size[1];
      }
    }

    empty_non_axis_sizes = (boffset == 0);
    if (empty_non_axis_sizes) {
      coffset = (int)(((1.0 + (double)i) - 1.0) * (double)r1->size[0]);
    } else {
      coffset = (int)(((1.0 + (double)i) - 1.0) * (double)r1->size[0]);
      loop_ub = r2->size[1];
      if (!((coffset == 0) || (loop_ub == 0))) {
        coffset = (int)(((1.0 + (double)i) - 1.0) * (double)r1->size[0]);
      } else {
        coffset = 0;
      }
    }

    if (empty_non_axis_sizes || (!((CAB->size[0] == 0) || (CAB->size[1] == 0))))
    {
      loop_ub = CAB->size[0];
    } else {
      loop_ub = 0;
    }

    for (i2 = 0; i2 < boffset; i2++) {
      for (c_i = 0; c_i < coffset; c_i++) {
        GAMMA->data[c_i + GAMMA->size[0] * ((i1 + i2) - 1)] = varargin_1->
          data[c_i + coffset * i2];
      }
    }

    for (i2 = 0; i2 < boffset; i2++) {
      for (c_i = 0; c_i < loop_ub; c_i++) {
        GAMMA->data[(c_i + coffset) + GAMMA->size[0] * ((i1 + i2) - 1)] =
          CAB->data[c_i + loop_ub * i2];
      }
    }

    i++;
  }

  emxFree_real_T(&varargin_1);
  emxFree_real_T(&CAB);
  emxFree_real_T(&y);
  emxFree_real_T(&r2);
  emxFree_real_T(&r1);
}

void mpc(const emxArray_real_T *A, const emxArray_real_T *B, const
         emxArray_real_T *C, const emxArray_real_T *x, double N, double r,
         double lb, emxArray_real_T *U)
{
  emxArray_real_T *PHI;
  int i0;
  int loop_ub;
  int i;
  emxArray_real_T *GAMMA;
  emxArray_real_T *r0;
  double u;
  int m;
  int aoffset;
  int inner;
  emxArray_real_T *a;
  emxArray_real_T *b;
  int coffset;
  int boffset;
  double ub;
  int b_i;
  int k;
  double y;
  emxInit_real_T(&PHI, 2);

  /*  Find matrix */
  /*  Create the special Observabillity matrix */
  i0 = PHI->size[0] * PHI->size[1];
  PHI->size[0] = (int)((double)C->size[0] * N);
  PHI->size[1] = A->size[1];
  emxEnsureCapacity_real_T(PHI, i0);
  loop_ub = (int)((double)C->size[0] * N) * A->size[1];
  for (i0 = 0; i0 < loop_ub; i0++) {
    PHI->data[i0] = 0.0;
  }

  i = 0;
  emxInit_real_T(&GAMMA, 2);
  emxInit_real_T(&r0, 2);
  while (i <= (int)N - 1) {
    u = (double)C->size[0] * (1.0 + (double)i);
    if (1.0 + (double)i > u) {
      i0 = 1;
    } else {
      i0 = i + 1;
    }

    mpower(A, 1.0 + (double)i, GAMMA);
    if ((C->size[1] == 1) || (GAMMA->size[0] == 1)) {
      aoffset = r0->size[0] * r0->size[1];
      r0->size[0] = C->size[0];
      r0->size[1] = GAMMA->size[1];
      emxEnsureCapacity_real_T(r0, aoffset);
      loop_ub = C->size[0];
      for (aoffset = 0; aoffset < loop_ub; aoffset++) {
        coffset = GAMMA->size[1];
        for (boffset = 0; boffset < coffset; boffset++) {
          r0->data[aoffset + r0->size[0] * boffset] = 0.0;
          b_i = C->size[1];
          for (k = 0; k < b_i; k++) {
            r0->data[aoffset + r0->size[0] * boffset] += C->data[aoffset +
              C->size[0] * k] * GAMMA->data[k + GAMMA->size[0] * boffset];
          }
        }
      }
    } else {
      m = C->size[0];
      inner = C->size[1];
      aoffset = r0->size[0] * r0->size[1];
      r0->size[0] = C->size[0];
      r0->size[1] = GAMMA->size[1];
      emxEnsureCapacity_real_T(r0, aoffset);
      for (loop_ub = 0; loop_ub < GAMMA->size[1]; loop_ub++) {
        coffset = loop_ub * m - 1;
        boffset = loop_ub * inner - 1;
        for (b_i = 1; b_i <= m; b_i++) {
          r0->data[coffset + b_i] = 0.0;
        }

        for (k = 1; k <= inner; k++) {
          if (GAMMA->data[boffset + k] != 0.0) {
            aoffset = (k - 1) * m;
            for (b_i = 1; b_i <= m; b_i++) {
              r0->data[coffset + b_i] += GAMMA->data[boffset + k] * C->data
                [(aoffset + b_i) - 1];
            }
          }
        }
      }
    }

    loop_ub = r0->size[1];
    for (aoffset = 0; aoffset < loop_ub; aoffset++) {
      coffset = r0->size[0];
      for (boffset = 0; boffset < coffset; boffset++) {
        PHI->data[((i0 + boffset) + PHI->size[0] * aoffset) - 1] = r0->
          data[boffset + r0->size[0] * aoffset];
      }
    }

    i++;
  }

  emxFree_real_T(&r0);
  gammaMat(A, B, C, N, GAMMA);

  /*  Solve first with no constraints */
  /*  Set U */
  i0 = U->size[0];
  U->size[0] = (int)N;
  emxEnsureCapacity_real_T1(U, i0);
  loop_ub = (int)N;
  for (i0 = 0; i0 < loop_ub; i0++) {
    U->data[i0] = 0.0;
  }

  /*  Iterate Gaussian Elimination */
  i = 0;
  emxInit_real_T(&a, 2);
  emxInit_real_T1(&b, 1);
  while (i <= (int)N - 1) {
    /*  Solve u */
    if (1.0 + (double)i == 1.0) {
      loop_ub = PHI->size[1];
      i0 = a->size[0] * a->size[1];
      a->size[0] = 1;
      a->size[1] = loop_ub;
      emxEnsureCapacity_real_T(a, i0);
      for (i0 = 0; i0 < loop_ub; i0++) {
        a->data[a->size[0] * i0] = PHI->data[PHI->size[0] * i0];
      }

      i0 = PHI->size[1];
      if ((i0 == 1) || (x->size[0] == 1)) {
        u = 0.0;
        for (i0 = 0; i0 < a->size[1]; i0++) {
          u += a->data[a->size[0] * i0] * x->data[i0];
        }
      } else {
        u = 0.0;
        for (i0 = 0; i0 < a->size[1]; i0++) {
          u += a->data[a->size[0] * i0] * x->data[i0];
        }
      }

      u = (r - u) / GAMMA->data[0];
    } else {
      loop_ub = PHI->size[1];
      i0 = a->size[0] * a->size[1];
      a->size[0] = 1;
      a->size[1] = loop_ub;
      emxEnsureCapacity_real_T(a, i0);
      for (i0 = 0; i0 < loop_ub; i0++) {
        a->data[a->size[0] * i0] = PHI->data[i + PHI->size[0] * i0];
      }

      i0 = PHI->size[1];
      if ((i0 == 1) || (x->size[0] == 1)) {
        u = 0.0;
        for (i0 = 0; i0 < a->size[1]; i0++) {
          u += a->data[a->size[0] * i0] * x->data[i0];
        }
      } else {
        u = 0.0;
        for (i0 = 0; i0 < a->size[1]; i0++) {
          u += a->data[a->size[0] * i0] * x->data[i0];
        }
      }

      loop_ub = (int)((1.0 + (double)i) - 1.0);
      i0 = a->size[0] * a->size[1];
      a->size[0] = 1;
      a->size[1] = (int)((1.0 + (double)i) - 1.0);
      emxEnsureCapacity_real_T(a, i0);
      for (i0 = 0; i0 < loop_ub; i0++) {
        a->data[a->size[0] * i0] = GAMMA->data[i + GAMMA->size[0] * i0];
      }

      loop_ub = (int)((1.0 + (double)i) - 1.0);
      i0 = b->size[0];
      b->size[0] = (int)((1.0 + (double)i) - 1.0);
      emxEnsureCapacity_real_T1(b, i0);
      for (i0 = 0; i0 < loop_ub; i0++) {
        b->data[i0] = U->data[i0];
      }

      if ((1.0 + (double)i) - 1.0 == 1.0) {
        y = 0.0;
        for (i0 = 0; i0 < a->size[1]; i0++) {
          y += a->data[a->size[0] * i0] * b->data[i0];
        }
      } else {
        y = 0.0;
        for (i0 = 0; i0 < a->size[1]; i0++) {
          y += a->data[a->size[0] * i0] * b->data[i0];
        }
      }

      u = ((r - u) - y) / GAMMA->data[i + GAMMA->size[0] * i];
    }

    /*  Constraints  */
    /*  Save u */
    U->data[i] = u;
    i++;
  }

  /*  Then use the last U as upper bound */
  ub = U->data[U->size[0] - 1];

  /*  Set U */
  i0 = U->size[0];
  U->size[0] = (int)N;
  emxEnsureCapacity_real_T1(U, i0);
  loop_ub = (int)N;
  for (i0 = 0; i0 < loop_ub; i0++) {
    U->data[i0] = 0.0;
  }

  /*  Iterate Gaussian Elimination */
  for (i = 0; i < (int)N; i++) {
    /*  Solve u */
    if (1.0 + (double)i == 1.0) {
      loop_ub = PHI->size[1];
      i0 = a->size[0] * a->size[1];
      a->size[0] = 1;
      a->size[1] = loop_ub;
      emxEnsureCapacity_real_T(a, i0);
      for (i0 = 0; i0 < loop_ub; i0++) {
        a->data[a->size[0] * i0] = PHI->data[PHI->size[0] * i0];
      }

      i0 = PHI->size[1];
      if ((i0 == 1) || (x->size[0] == 1)) {
        u = 0.0;
        for (i0 = 0; i0 < a->size[1]; i0++) {
          u += a->data[a->size[0] * i0] * x->data[i0];
        }
      } else {
        u = 0.0;
        for (i0 = 0; i0 < a->size[1]; i0++) {
          u += a->data[a->size[0] * i0] * x->data[i0];
        }
      }

      u = (r - u) / GAMMA->data[0];
    } else {
      loop_ub = PHI->size[1];
      i0 = a->size[0] * a->size[1];
      a->size[0] = 1;
      a->size[1] = loop_ub;
      emxEnsureCapacity_real_T(a, i0);
      for (i0 = 0; i0 < loop_ub; i0++) {
        a->data[a->size[0] * i0] = PHI->data[i + PHI->size[0] * i0];
      }

      i0 = PHI->size[1];
      if ((i0 == 1) || (x->size[0] == 1)) {
        u = 0.0;
        for (i0 = 0; i0 < a->size[1]; i0++) {
          u += a->data[a->size[0] * i0] * x->data[i0];
        }
      } else {
        u = 0.0;
        for (i0 = 0; i0 < a->size[1]; i0++) {
          u += a->data[a->size[0] * i0] * x->data[i0];
        }
      }

      loop_ub = (int)((1.0 + (double)i) - 1.0);
      i0 = a->size[0] * a->size[1];
      a->size[0] = 1;
      a->size[1] = (int)((1.0 + (double)i) - 1.0);
      emxEnsureCapacity_real_T(a, i0);
      for (i0 = 0; i0 < loop_ub; i0++) {
        a->data[a->size[0] * i0] = GAMMA->data[i + GAMMA->size[0] * i0];
      }

      loop_ub = (int)((1.0 + (double)i) - 1.0);
      i0 = b->size[0];
      b->size[0] = (int)((1.0 + (double)i) - 1.0);
      emxEnsureCapacity_real_T1(b, i0);
      for (i0 = 0; i0 < loop_ub; i0++) {
        b->data[i0] = U->data[i0];
      }

      if ((1.0 + (double)i) - 1.0 == 1.0) {
        y = 0.0;
        for (i0 = 0; i0 < a->size[1]; i0++) {
          y += a->data[a->size[0] * i0] * b->data[i0];
        }
      } else {
        y = 0.0;
        for (i0 = 0; i0 < a->size[1]; i0++) {
          y += a->data[a->size[0] * i0] * b->data[i0];
        }
      }

      u = ((r - u) - y) / GAMMA->data[i + GAMMA->size[0] * i];
    }

    /*  Constraints  */
    if (u > ub) {
      u = ub;
    } else {
      if (u < lb) {
        u = lb;
      }
    }

    /*  Save u */
    U->data[i] = u;
  }

  emxFree_real_T(&b);
  emxFree_real_T(&a);
  emxFree_real_T(&GAMMA);
  emxFree_real_T(&PHI);
}

/* End of code generation (mpc.c) */

MATLAB filen är ca 70 rader kod med lite for-loopar och if-satser.

[Icecap]

Jag har aldrig använd dynamisk allokering under alla de år jag har programmerat µC!

Jag har däremot ofta allokerat buffrar i programmen och inte sällan återanvänd dom (efter NOGA analys!) vid att deklarera dom som union{}.

Men efter en snabb titt ser det ut att det program du visar fram använder recursive programmering och då kan det vara nödvändigt kan jag tänka mig.

Jag antar att det är en µC med en del resurser och så kan man ju ibland tänka på det hela lite som ett SoC - där det kan vara vettigt att ha dynamisk allokering.

Men om en manuellt genererat kod blir signifikant mindre betyder det ofta också att den kör en del snabbare - något som man ska väga med i det hela.

[TomasL]

Dynamisk allokering är väl inte tillåtet enligt MISRA och andra standarder.
Ej heller rekursiva funktioner.

[DanielM]

Jag har aldrig använd dynamisk allokering under alla de år jag har programmerat µC!

Jag har däremot ofta allokerat buffrar i programmen och inte sällan återanvänd dom (efter NOGA analys!) vid att deklarera dom som union{}.

Men efter en snabb titt ser det ut att det program du visar fram använder recursive programmering och då kan det vara nödvändigt kan jag tänka mig.

Jag antar att det är en µC med en del resurser och så kan man ju ibland tänka på det hela lite som ett SoC - där det kan vara vettigt att ha dynamisk allokering.

Men om en manuellt genererat kod blir signifikant mindre betyder det ofta också att den kör en del snabbare - något som man ska väga med i det hela.

Jag postar koden här.

Kod: Markera allt

function [U] = mpc (A, B, C, x, N, r, lb)

  % Find matrix
  PHI = phiMat(A, C, N)
  GAMMA = gammaMat(A, B, C, N);

  % Solve first with no constraints
  U = solve(PHI, GAMMA, x, N, r, 0, 0, false);
  
  % Then use the last U as upper bound
  U = solve(PHI, GAMMA, x, N, r, lb, U(end), true);
  
end

function U = solve(PHI, GAMMA, x, N, r, lb, ub, constraints)
  % Set U
  U = zeros(N, 1);
  
  % Iterate Gaussian Elimination
  for i = 1:N
    
    % Solve u
    if(i == 1)
      u = (r - PHI(i,:)*x)/GAMMA(i,i);
    else
      u = (r - PHI(i,:)*x - GAMMA(i,1:i-1)*U(1:i-1) )/GAMMA(i,i);
    end
    
    % Constraints 
    if(constraints == true)
      if(u > ub)
        u = ub;
      elseif(u < lb)
        u = lb;
      end
    end
    
    % Save u
    U(i) = u
  end
end

function PHI = phiMat(A, C, N)

  % Create the special Observabillity matrix
  PHI = zeros(size(C,1)*N, size(A, 2));
  for i = 1:N
    PHI(i:size(C,1)*i, :) = C*A^i;
  end

end

function GAMMA = gammaMat(A, B, C, N)

  % Create the lower triangular toeplitz matrix
  GAMMA = zeros(size(C,1)*N, size(B, 2)*N);
  for i = 1:N
    GAMMA(:, i:size(C,1)*i) = vertcat(zeros((i-1)*size(C*A*B, 1), size(C*A*B, 2)),cabMat(A, B, C, N-i+1));
  end

end

function CAB = cabMat(A, B, C, N)

  % Create the column for the GAMMA matrix
  CAB = zeros(size(C,1)*N, size(B, 2));
  for i = 0:N-1
    CAB((i+1):size(C,1)*(i+1), :) = C*A^i*B;
  end

end

Som du ser så är den svåra matematiken följande:

Kod: Markera allt

A^i
vertcat

Dynamisk allokering är väl inte tillåtet enligt MISRA och andra standarder.
Ej heller rekursiva funktioner.

Varför inte då?
Detta är ett reglersystem som jag håller på med.

[TomasL]

Av den enkla anledningen att det blir buggar.
Man behöver varken rekursiva funktioner eller dynamisk allokering för ett reglersystem.

[DanielM]

Jag köper att koden blir komplex. Men borde inte MATLAB's C-coder utvecklingsteam vara medveten om detta? MATLAB har väldigt många toolboxar för att generera C-kod för inbyggda system.

Tänkte då passar jag på och checka läget med deras MATLAB Coder.

Edit:
Såg nu att MATLAB Coder har val att inte generera dynamiskt minne.

[Gimbal]

När jag jobbade med inbyggda system så var minnesallokeringar endast tillåtet under uppstart, sedan fick det inte röras. Man ville inte riskera minnesfragmentering för allt smör i småland.

[Krille Krokodil]

Numerisk analys utan minnesallokering kan vara en liten huvudvärk att få till men det skulle kanske gå att ta
någon mellanväg där man gör egen mindre generell minneshantering där man har förallokerat minne till ett
max antal matriser samt har en felhantering om Matlab-koden skulle begära fler än den övre gräns man har
satt.

[DanielM]

När jag jobbade med inbyggda system så var minnesallokeringar endast tillåtet under uppstart, sedan fick det inte röras. Man ville inte riskera minnesfragmentering för allt smör i småland.

Kanske man ska tillämpa structs i MATLAB för att ha matriserna där i?

Numerisk analys utan minnesallokering kan vara en liten huvudvärk att få till men det skulle kanske gå att ta
någon mellanväg där man gör egen mindre generell minneshantering där man har förallokerat minne till ett
max antal matriser samt har en felhantering om Matlab-koden skulle begära fler än den övre gräns man har
satt.

Mycket svårt kan jag säga! Jag har försökt skriva egen C-kod, men det är många timmars spenderande för att få beräkningarna rätt.

[DanielM]

Efter att ha spenderat några timmar så fungerar denna kod exakt som MATLAB filen ovan. Ge gärna kommentarer på den. Jag vet att det går alltid göra den bättre, jag har själv några förslag också för att snygga till den. Men fler ögon är bättre än två. els har jag inte jobbat klart på än.

Jag har som mål att denna ska passa alla inbyggda system, oavsett CPU eller RAM. Jag ska minimiera så mycket som det bara går och göra det enkelt.

Kod: Markera allt

/*
 * Model_Predictive_Control.c
 *
 *  Created on:
 *      Author:
 */

#include "Model_Predictive_Control.h"

/*
 * Parameters
 */
int adim;
int ydim;
int rdim;
int horizon;

/*
 * Deceleration
 */
static void obsv(float* PHI, float* A, float* C);
static void kalman(float* x, float* A, float* B, float* u, float* K, float* y, float* C);
static void mul(float* A, float* B, float* C, int row_a, int column_a, int column_b);
static void tran(float* A, int row, int column);
static void CAB(float* GAMMA, float* PHI, float* A, float* B, float* C);
static void solve(float* GAMMA, float* PHI, float* x, float* u, float* r, float lb, float ub, int constraintsON);
static void print(float* A, int row, int column);


void mpc(int adim_, int ydim_, int rdim_, int horizon_, float* A, float* B, float* C, float* D, float* K, float* u, float* r, float* y, float* x){
	/*
	 * Set the dimensions
	 */
	adim = adim_;
	ydim = ydim_;
	rdim = rdim_;
	horizon = horizon_;

	/*
	 * Identify the model - Extended Least Square
	 */
	int n = 5;
	float* phi;
	float* theta;
	//els(phi, theta, n, y, u, P);

	/*
	 * Create a state space model with Observable canonical form
	 */


	/*
	 * Create the extended observability matrix
	 */
	float PHI[horizon*ydim*adim];
	memset(PHI, 0, horizon*ydim*adim*sizeof(float));
	obsv(PHI, A, C);

	/*
	 * Create the lower triangular toeplitz matrix
	 */
	float GAMMA[horizon*rdim*horizon*ydim];
	memset(GAMMA, 0, horizon*rdim*horizon*ydim*sizeof(float));
	CAB(GAMMA, PHI, A, B, C);

	/*
	 * Solve the best input value
	 */
	solve(GAMMA, PHI, x, u, r, 0, 0, 0);
	solve(GAMMA, PHI, x, u, r, 0, *(u), 1);

	/*
	 * Estimate the state vector
	 */
	kalman(x, A, B, u, K, y, C);
}

/*
 * Identify the model
 */
static void els(float* P, float* phi, float* theta, int polyLength, int totalPolyLength, float* y, float* u, float* e){

	/*
	 * move phi with the inputs, outputs, errors one step to right
	 */
	for(int i = 0; i < polyLength; i++){
		*(phi + i+1 + totalPolyLength*0) = *(phi + i + totalPolyLength*0); // Move one to right for the y's
		*(phi + i+1 + totalPolyLength*1) = *(phi + i + totalPolyLength*1); // Move one to right for the u's
		*(phi + i+1 + totalPolyLength*2) = *(phi + i + totalPolyLength*2); // Move one to right for the e's
	}

	/*
	 * Add the current y, u and e

	(*phi + totalPolyLength*0) = -*(y + 0); // Need to be negative!
	(*phi + totalPolyLength*1) =  *(u + 0);
	(*phi + totalPolyLength*2) =  *(e + 0);
 */
	/*
	 * phi'*theta
	 */
	float y_est = 0;
	for(int i = 0; i < totalPolyLength; i++){
		y_est += *(phi + i) * *(theta + i);
	}
	float epsilon = *(y + 0) - y_est; // In this case, y is only one element array

	/*
	 * phi*epsilon
	 */
	float phi_epsilon[totalPolyLength];
	memset(phi_epsilon, 0, totalPolyLength*sizeof(float));
	for(int i = 0; i < totalPolyLength; i++){
		*(phi_epsilon + i) = *(phi + i) * epsilon;
	}

	/*
	 * P_vec = P*phi_epsilon
	 */
	float P_vec[totalPolyLength];
	memset(P_vec, 0, totalPolyLength*sizeof(float));
	mul(P, phi_epsilon, P_vec, totalPolyLength, totalPolyLength, 1);

	/*
	 * Update our estimated vector theta = theta + P_vec
	 */
	for(int i = 0; i < totalPolyLength; i++){
		*(theta + i) = *(theta + i) + *(P_vec + i);
	}

	/*
	 * Update P = P - (P*phi*phi'*P)/(1 + phi'*P*phi)
	 */
	// Create phi'
	float phiT[totalPolyLength];
	memset(phiT, 0, totalPolyLength*sizeof(float));
	memcpy(phiT, phi, totalPolyLength*sizeof(float));
	tran(phiT, totalPolyLength, 1);

	// phi'*P
	float phiT_P[totalPolyLength];
	memset(phiT_P, 0, totalPolyLength*sizeof(float));
	mul(phiT, P, phiT_P, 1, totalPolyLength, totalPolyLength);

	// phi*phi'*P
	float phi_phiT_P[totalPolyLength*totalPolyLength];
	memset(phi_phiT_P, 0, totalPolyLength*totalPolyLength*sizeof(float));
	mul(phi, phiT_P, phi_phiT_P, totalPolyLength, 1, totalPolyLength);

	// P*phi*phi'*P
	float P_phi_phiT_P[totalPolyLength*totalPolyLength];
	memset(P_phi_phiT_P, 0, totalPolyLength*totalPolyLength*sizeof(float));
	mul(P, phi_phiT_P, P_phi_phiT_P, totalPolyLength, totalPolyLength, totalPolyLength);

	// P*phi
	float P_phi[totalPolyLength];
	memset(P_phi, 0, totalPolyLength*sizeof(float));
	mul(P, phi, P_phi, totalPolyLength, totalPolyLength, 1);

	// phi'*P*phi
	float phiT_P_phi[1];
	memset(phiT_P_phi, 0, 1*sizeof(float));
	mul(phiT, P_phi, phiT_P_phi, 1, totalPolyLength, 1);

	// P = P - (P_phi_phiT_P) / (1+phi'*P*phi)
	for(int i = 0; i < totalPolyLength*totalPolyLength; i++){
		*(P + i) = *(P + i) - *(P_phi_phiT_P + i) / (1 + *(phiT_P_phi));
	}

}

/*
 * This will solve if GAMMA is square!
 */
static void solve(float* GAMMA, float* PHI, float* x, float* u, float* r, float lb, float ub, int constraintsON){
	/*
	 * Now we are going to solve on the form
	 * Ax=b, where b = (R*r-PHI*x) and A = GAMMA and x = U
	 */

	/*
	 * R_vec = R*r
	 */
	float R_vec[horizon*ydim];
	memset(R_vec, 0, horizon*ydim*sizeof(float));
	for(int i = 0; i < horizon*ydim; i++){
		for (int j = 0; j < rdim; j++) {
			*(R_vec + i + j) = *(r + j);
		}
		i += rdim-1;
	}

	/*
	 * PHI_vec = PHI*x
	 */
	float PHI_vec[horizon*ydim];
	memset(PHI_vec, 0, horizon * ydim * sizeof(float));
	mul(PHI, x, PHI_vec, horizon*ydim, adim, 1);

	/*
	 * Solve now (R_vec - PHI_vec) = GAMMA*U
	 * Notice that this is ONLY for Square GAMMA with lower triangular toeplitz matrix e.g SISO case
	 * This using Gaussian Elimination backward substitution
	 */
	float U[horizon];
	float sum = 0.0;
	memset(U, 0, horizon*sizeof(float));
	for(int i = 0; i < horizon; i++){
		for(int j = 0; j < i; j++){
			sum += *(GAMMA + i*horizon + j) * *(U + j);
		}
		float newU = (*(R_vec + i) - *(PHI_vec + i) - sum) / (*(GAMMA + i*horizon + i));
		if(constraintsON == 1){
			if(newU > ub)
				newU = ub;
			if(newU < lb)
				newU = lb;
		}
		*(U + i) = newU;
		sum = 0.0;
	}
	//print(U, horizon, 1);

	/*
	 * Set last U to u
	 */
	if(constraintsON == 0){
		*(u + 0) = *(U + horizon - 1);
	}else{
		*(u + 0) = *(U + 0);
	}


}

/*
 * Lower traingular toeplitz of extended observability matrix
 */
static void CAB(float* GAMMA, float* PHI, float* A, float* B, float* C){
	/*
	 * First create the initial C*A^0*B == C*I*B == C*B
	 */
	float CB[ydim*rdim];
	memset(CB, 0, ydim*rdim*sizeof(float));
	mul(C, B, CB, ydim, adim, rdim);

	/*
	 * Take the transpose of CB so it will have dimension rdim*ydim instead
	 */
	tran(CB, ydim, rdim);

	/*
	 * Create the CAB matrix from PHI*B
	 */
	float PHIB[horizon*ydim*rdim];
	mul(PHI, B, PHIB, horizon*ydim, adim, rdim); // CAB = PHI*B
	tran(PHIB, horizon*ydim, rdim);

	/*
	 * We insert GAMMA = [CB PHI;
	 *                    0  CB PHI;
	 *            		  0   0  CB PHI;
	 *            		  0   0   0  CB PHI] from left to right
	 */
	for(int i = 0; i < horizon; i++) {
		for(int j = 0; j < rdim; j++) {
			memcpy(GAMMA + horizon*ydim*(i*rdim+j) + ydim*i, CB + ydim*j, ydim*sizeof(float)); // Add CB
			memcpy(GAMMA + horizon*ydim*(i*rdim+j) + ydim*i + ydim, PHIB + horizon*ydim*j, (horizon-i-1)*ydim*sizeof(float)); // Add PHI*B
		}
	}
	/*
	 * Transpose of gamma
	 */
	tran(GAMMA, horizon*rdim, horizon*ydim);
	//print(CB, rdim, ydim);
	//print(PHIB, rdim, horizon*ydim);
	//print(GAMMA, horizon*ydim, horizon*rdim);
}

/*
 * Transpose
 */
static void tran(float* A, int row, int column) {

	float B[row*column];
	float* transpose;
	float* ptr_A = A;

	for (int i = 0; i < row; i++) {
		transpose = &B[i];
		for (int j = 0; j < column; j++) {
			*transpose = *ptr_A;
			ptr_A++;
			transpose += row;
		}
	}

	// Copy!
	memcpy(A, B, row*column*sizeof(float));

}

/*
 * [C*A^1; C*A^2; C*A^3; ... ; C*A^horizon] % Extended observability matrix
 */
static void obsv(float* PHI, float* A, float* C){

	/*
	 * This matrix will A^(i+1) all the time
	 */
	float A_pow[adim*adim];
	memset(A_pow, 0, adim * adim * sizeof(float));
	float A_copy[adim*adim];
	memcpy(A_copy, A, adim * adim * sizeof(float));

	/*
	 * Temporary matrix
	 */
	float T[ydim*adim];
	memset(T, 0, ydim * adim * sizeof(float));

	/*
	 * Regular T = C*A^(1+i)
	 */
	mul(C, (float*) A, T, ydim, adim, adim);

	/*
	 * Insert temporary T into PHI
	 */
	memcpy(PHI, T, ydim*adim*sizeof(float));

	/*
	 * Do the rest C*A^(i+1) because we have already done i = 0
	 */
	for(int i = 1; i < horizon; i++){
		mul(A, A_copy, A_pow, adim, adim, adim); //  Matrix power A_pow = A*A_copy
		mul(C, A_pow, T, ydim, adim, adim); // T = C*A^(1+i)
		memcpy(PHI + i*ydim*adim, T, ydim*adim*sizeof(float)); // Insert temporary T into PHI
		memcpy(A_copy, A_pow, adim * adim * sizeof(float)); // A_copy <- A_pow
	}
}

/*
 *   x = Ax - KCx + Bu + Ky % Kalman filter
 */
static void kalman(float* x, float* A, float* B, float* u, float* K, float* y, float* C) {
	/*
	 * Compute the vector A_vec = A*x
	 */
	float A_vec[adim*1];
	memset(A_vec, 0, adim*sizeof(float));
	mul(A, x, A_vec, adim, adim, 1);

	/*
	 * Compute the vector B_vec = B*u
	 */
	float B_vec[adim*1];
	memset(B_vec, 0, adim*sizeof(float));
	mul(B, u, B_vec, adim, rdim, 1);

	/*
	 * Compute the vector C_vec = C*x
	 */
	float C_vec[ydim*1];
	memset(C_vec, 0, ydim*sizeof(float));
	mul(C, x, C_vec, ydim, adim, 1);

	/*
	 * Compute the vector KC_vec = K*C_vec
	 */
	float KC_vec[adim*1];
	memset(KC_vec, 0, adim*sizeof(float));
	mul(K, C_vec, KC_vec, adim, ydim, 1);

	/*
	 * Compute the vector Ky_vec = K*y
	 */
	float Ky_vec[adim*1];
	memset(Ky_vec, 0, adim*sizeof(float));
	mul(K, y, Ky_vec, adim, ydim, 1);

	/*
	 * Now add x = A_vec - KC_vec + B_vec + Ky_vec
	 */
	for(int i = 0; i < adim; i++){
		*(x + i) = *(A_vec + i) - *(KC_vec + i) + *(B_vec + i) + *(Ky_vec + i);
	}

}

/*
 * C = A*B
 */
static void mul(float* A, float* B, float* C, int row_a, int column_a, int column_b) {

	// Data matrix
	float* data_a = A;
	float* data_b = B;

	for (int i = 0; i < row_a; i++) {
		// Then we go through every column of b
		for (int j = 0; j < column_b; j++) {
			data_a = &A[i * column_a];
			data_b = &B[j];

			*C = 0; // Reset
			// And we multiply rows from a with columns of b
			for (int k = 0; k < column_a; k++) {
				*C += *data_a * *data_b;
				data_a++;
				data_b += column_b;
			}
			C++; // ;)
		}
	}
}

/*
 * Print matrix or vector - Just for error check
 */
static void print(float* A, int row, int column) {
	for (int i = 0; i < row; i++) {
		for (int j = 0; j < column; j++) {
			printf("%0.18f ", *(A++));
		}
		printf("\n");
	}
	printf("\n");
}

Med GCC: 0.000087 sekunder
Med Octave openBLAS: 0.0362771 sekunder

Om ni undrar vad jag sjutton håller på med så kretsar mina dagar och nätter kring detta. Detta är en tillståndsmodell, en linjär sådan.
\(x(k+1) = Ax(k) + Bu(k)\)
\(y(k) = Cx(k)\)

Fördelen men modellen är att vi kan prediktivt "förutspå" framtiden lite enkelt förklarat.

\(\Phi = \begin{bmatrix}
CA\\
CA^2\\
CA^3\\
\vdots \\
CA^{n-1}
\end{bmatrix}\)

\(\Gamma = \begin{bmatrix}
CB & 0 & 0 & 0 &0 \\
CAB & CB & 0 & 0 &0 \\
CA^2B & CAB & CB & 0 & 0\\
\vdots & \vdots &\vdots & \ddots & 0 \\
CA^{n-2} B& CA^{n-3} B & CA^{n-4} B & \dots & CA^{n-j} B
\end{bmatrix}\)

Där \(j=n\)

Då löser vi detta igenom minstakvadratmetoden. Vi vill alltså veta "Vad ska jag ha för insignal till systemet för att systemet ska vara prediktiv?".

\(U = (\Gamma)^{-1}(R-\Phi x)\)

Men tittar vi på insignalen \(U\) så beter den sig så här.


Den blåa signalen är dead-beat reglering. Ser ut så här i en PC simulering.


Men använder vi min idé ovan så får vi ett signalsystem som ser ut så här:


Vilket resulterar ett mer verkligt beteende.

[DanielM]

Jag uppskattar gärna om jag skulle kunna få några tips&tricks angående någon hård C-standard så som MIRSA C. Vad ska jag tänka på om jag vill skriva säker C-kod?

[hummel]

MISRA C är en bra standard, det finns förövrigt 3 versioner. Jag har själv varit ansvarig för införandet av MISRA C:1998 på ett företag för fordonselektronik.
Det knepiga är att förstå tänket bakom, varför finns dessa regler och vad fyller dom för funktion. En regel säger att assembler-kod inte är tillåtet. I vissa implementeringar kan det var god kodning att tillåta ett par rader assembler, det är val man själv får göra i projektet. Du kan göra ett projekt som följer MISRA C:1998 med undantag för vissa av dess 127 regler.
Det är tänket som är viktigast (som alltid). Dom flesta reglerna i MISRA C är rätt självklara om man arbetat med inbyggda system med säkerhetskritisk funktionalitet. Fördelen är att man slipper komma på allt själv.

PC-Lint från Gimpel är ett suveränt verktyg för statisk analys av C-kod (finner buggar man själv nätt och jämt lyckats skapa. ) som kan konfigureras för MISRA C.

[DanielM]

Trodde MIRSA C var ett regelverk som endast förbjöd dynamiskt allokering utav minne. Visste inte att Assembler var en säkerhetsrisk.

[sodjan]

Inte senaste versionen, men det ger väl en klar indikering på
att det inte enbart gäller allokering av minne...

http://caxapa.ru/thumbs/468328/misra-c-2004.pdf

[DanielM]

Antar att du är väl bekant med denna standard? Skulle du kunna peka på operationer i min C kod där jag gör fel enligt standarden?