复习了一下优化的主要形式, 这里关于一些重点做一下笔记.

关于二次规划的H矩阵

二次规划的标准形式如下.

在matlab中二次规划可以用quadprog, 一次规划用linprog, 0-1规划用bintprog, 最简单的就是fminbnd了, 另外常用的就是fminunc, fmincon, fzero, 对最大最小问题, 一般转换成最小问题解, 这在mathematica和lingo里面都是这样处理的, 但是matlab里面有内置的函数fminmax, 第一个参数为一个行向量, 为取最大的一个, 多目标用fgoalattain, 这个感觉实用性不是很大.
对于quadprog, matlab里面调用方法如下.

主要是这里标准形式的H矩阵有点难理解.
以一个例子来看

其解法如下

其中的f向量比较好理解, 对于H矩阵, 其具体求法去下.

所以剩下的简单了.

题目

主要题意就是在地下水管网络上建设蓄水池来缓解地下水道的排水压力, 如下图.

储水池的参数如下表所示，试对该旧城区地下管网改建工程问题进行建模，并提供一种可行的改造方案。

上表就是题目给的唯一数据….什么鬼, 给个地图不给海拔经纬度信息就给个图是神马, 连个比例尺都没有..摔…

模型就不讨论了, 主要就是个优化模型问题,(上去就是一个蛙叫算法)哈哈..

黑魔法阶段

由于没有数据, 我就开始了搜索之旅…..结果居然用google map找到了题目中的区域…在葡萄牙的克英布拉, 哈哈…
google map上的该地点.

利用google map可以将这块cathchment围起来, 计算其总的面积和周长.

之后继续google altitude map, 可以把这17个storage unit的海拔取出来, 如下表

这题是一个光伏角的优化问题, 模型问题不再赘述, 我是负责算法部分, 所以在这里重述一下算法的思路.

题目

太阳能是一种清洁能源，它的应用正在世界范围内快速地增长。目前建设一个太阳能发电系统成本较高，从我国现阶段的太阳能发电成本来看，其花费在太阳电池组件的费用大约为60％左右。为了更好地利用太阳能电池板发电，从安装的角度来看，也有可以挖掘的潜力。

涉及到的几何参数

1. 太阳能电池板的倾角（Module tilt）：太阳能电池板平面与水平面的夹角
2. 太阳能电池板的方位角（Module azimuth）：此角度为太阳能电池板的法线在水平面上的投影与正南方向的夹角。并按如下方法规定角度：如果一个太阳能电池板朝向正南，那么它的方位角为零度; 如果一个太阳能电池板朝向正东，那么它的方位角为-90度; 如果一个太阳能电池板朝向正西，那么它的方位角为+90度。
3. 太阳高度角：对于地球表面上的某个点，太阳高度角是指太阳光的入射方向和地平面之间的夹角。
4. 太阳方位角：在地球表面某点，建立以球面外法线方向为z轴正方向，以向东方向为x轴，向北方向为y轴的坐标系。太阳在这个坐标系中，北偏东的角度为方位角。
5. 经度、纬度、海拔高度

latitude： in degrees, north of equator is positive
longitude ：in degrees, positive for east of Greenwich
altitude: in meters ，above mean sea level

污水排放问题

这次出的题目是一个污水扩散的一个模型, 需要计算在下游3km处的一个工厂的排放污染量最大排多少才能使得在下游10km处的一个监测站监测的数据不超过10mg/L.

这次的模型相对复杂, 最后我们结果也不是太满意, 做了一个2维扩散的模型, 我是编程的, 对模型部分就不多做讨论了, 主要进行两个方面的讨论, 工厂是连续投放污染物还是一次性投放污染物, 对于一次投放的计算代码如下.

关于隐函数相交线问题我们数模队友问下了, 话说他每天弄他们校创飞机的事情搞笑的= =.进入正题, 关于两函数相交线, 我们只需要找到同时满足两函数方程的xyz即可

matlab中对于3d绘图有isosurface还有alphashape, boundary, trisurf, convhull, surf等, 每个函数各自实现的方法和效果都不太一样, 这里用isosurface进行一个说明, 其传入参数需要为一个3维数组, 绘制的结果为其中为0的值, 也就是满足函数关系等于0的值, 这里举个例子说明

关于 R 中 scoping rules 的一些笔记与理解.

在一个Function中，对一个变量的lexical scoping rule的检索次序是，先检查本Function中的environment，如果能找到该变量，则返回该变量，如果不能，则检索这个Function被创建(不是被调用)的Environment，以此类推，一直到检索Global Environment。

作为一个小总结吧…
学习是从上个寒假开始的, 熊老师让我帮他弄那个轮椅控制的事情= =当时每个星期天都跑了过去帮忙, 连续弄了好几个星期吧, 最后快到期末基本弄成了.

然后这两天帮他们把他们后续的东西大致给弄完了, 然后晚上感觉在那边睡着不是很踏实, 就今天跑回来了, 毕竟后面的任务我在家也可以完成了, 在那边确实呆着不是很舒服, 洗澡什么的也不方便.

帮那个同学弄了一下震荡检测的, 虽然这个算法我觉得还是需要在考虑下或者参数需要调好点, 不过最后结果看起来还是可以的.

在那边主要帮忙了三个东西, 一个摇晃检测, 一个移动动态控制, 还有一个GPRS数据传输的东西.

不多说代码如下

摇晃检测

#include <CurieIMU.h>
#include <MadgwickAHRS.h>

Madgwick filter; // initialise Madgwick object
int ax, ay, az;
int gx, gy, gz;
float yaw[30];
float pitch[30];
float roll[30];
int factor = 800; // variable by which to divide gyroscope values, used to control sensitivity
// note that an increased baud rate requires an increase in value of factor

int calibrateOffsets = 1; // int to determine whether calibration takes place or not


void setup() {
  // initialize Serial communication
  Serial.begin(9600);

  // initialize device
  CurieIMU.begin();

  if (calibrateOffsets == 1) {
    // use the code below to calibrate accel/gyro offset values
    Serial.println("Internal sensor offsets BEFORE calibration...");
    Serial.print(CurieIMU.getAccelerometerOffset(X_AXIS)); Serial.print("\t");
    Serial.print(CurieIMU.getAccelerometerOffset(Y_AXIS)); Serial.print("\t");
    Serial.print(CurieIMU.getAccelerometerOffset(Z_AXIS)); Serial.print("\t");
    Serial.print(CurieIMU.getGyroOffset(X_AXIS)); Serial.print("\t");
    Serial.print(CurieIMU.getGyroOffset(Y_AXIS)); Serial.print("\t");
    Serial.print(CurieIMU.getGyroOffset(Z_AXIS)); Serial.print("\t");
    Serial.println("");

    // To manually configure offset compensation values, use the following methods instead of the autoCalibrate...() methods below
    //    CurieIMU.setGyroOffset(X_AXIS, 220);
    //    CurieIMU.setGyroOffset(Y_AXIS, 76);
    //    CurieIMU.setGyroOffset(Z_AXIS, -85);
    //    CurieIMU.setAccelerometerOffset(X_AXIS, -76);
    //    CurieIMU.setAccelerometerOffset(Y_AXIS, -235);
    //    CurieIMU.setAccelerometerOffset(Z_AXIS, 168);

    //IMU device must be resting in a horizontal position for the following calibration procedure to work correctly!

    Serial.print("Starting Gyroscope calibration...");
    CurieIMU.autoCalibrateGyroOffset();
    Serial.println(" Done");
    Serial.print("Starting Acceleration calibration...");
    CurieIMU.autoCalibrateAccelerometerOffset(X_AXIS, 0);
    CurieIMU.autoCalibrateAccelerometerOffset(Y_AXIS, 0);
    CurieIMU.autoCalibrateAccelerometerOffset(Z_AXIS, 1);
    Serial.println(" Done");

    Serial.println("Internal sensor offsets AFTER calibration...");
    Serial.print(CurieIMU.getAccelerometerOffset(X_AXIS)); Serial.print("\t");
    Serial.print(CurieIMU.getAccelerometerOffset(Y_AXIS)); Serial.print("\t");
    Serial.print(CurieIMU.getAccelerometerOffset(Z_AXIS)); Serial.print("\t");
    Serial.print(CurieIMU.getAccelerometerOffset(X_AXIS)); Serial.print("\t");
    Serial.print(CurieIMU.getAccelerometerOffset(Y_AXIS)); Serial.print("\t");
    Serial.print(CurieIMU.getAccelerometerOffset(Z_AXIS)); Serial.print("\t");
    Serial.println("");
  }
}

int k = 0;
int flag = 0;
float roll3[3];
float pitch3[3];
int sign = 0;

void readFirstMotion(){
  for(int i = 0;i<30;i++){
    CurieIMU.readMotionSensor(ax, ay, az, gx, gy, gz); 
    filter.updateIMU(gx / factor, gy / factor, gz / factor, ax, ay, az);
    roll[i] = filter.getRoll();
    pitch[i] = filter.getPitch();
  }
}

void refreshMotion(){
  for(int i = 0;i<29;i++){
    roll[i] = roll[i+1];
    pitch[i] = pitch[i+1];
  }
  CurieIMU.readMotionSensor(ax, ay, az, gx, gy, gz); 
  filter.updateIMU(gx / factor, gy / factor, gz / factor, ax, ay, az);
  roll[29] = filter.getRoll();
  pitch[29] = filter.getPitch();
}


int var(){
  refreshMotion();
  for(int i = 0;i<29;i++){
    if(i % 10 == 0 && i <= 20){  
      float roll_mean = 0;
      float pitch_mean = 0;
      for(int k = i;k<i+10;k++){
        roll_mean += roll[k];
        pitch_mean += pitch[k];
      }
      int m = i/10;
      roll3[m] = roll_mean/10;
      pitch3[m] = pitch_mean/10;
    }
  }
 if(abs((roll3[1] - roll3[2]) >= 0.5) || (abs(roll3[2] - roll[3]) >= 0.5) || abs((pitch3[1] - pitch3[2]) >= 0.5) || (abs(pitch3[2] - pitch[3]) >= 0.5))
   return 1;
 else 
   return 0;
}

int judge(){
  int count = 0;
  for(int i = 0; i<5000; i++){
    int sign = var();
    if(sign == 1){
      count++;
    }
  }
  if(count >= 4000){
    return 1;
  }else
    return 0;
}

void loop() {
  if(flag == 0){
    readFirstMotion();
    flag++;
  }
  sign = judge();
  if (sign == 0){
    Serial.println("situation is OK");
  }else{
    Serial.println("bad situation");
  }
}

龙门吊问题

这次出的题是一个龙门吊运动过程的角度优化问题, 主要是多目标规划和微分方程, 感觉这个多目标规划的目标函数很难找到一个满意的, 几个目标的无量纲很难找到一个客观的评判标准, 一些参数不能自己主观设置, 但是又没有客观的方法, 确实很麻烦.

然后貌似是陈建业从网上找来的, 直接找到了原题= =不过没有结果, 然而我google搜了一下gantry crane….居然找到了mathworks的一个mathematical modeling的包….里面就包含了这个模型的仿真结果, 不过觉得第一次训练还是好好写写主要不是结果还是要磨合一下队伍就没有用那个结果, 也没和队友说. 自己重新做了一个..下面是用matlab做的一个结果, 感觉还行.

首先是一个随机生成的结果

这里用的是随机优化, 就是个random search….用的是matlab的fmincon..感觉这种题基本很难找到全局最优解, 每一步都有一个微分方程要解, 估计不可导, 就是这个原因弄的我第一天弄了一下午的时间改mathematica的NMinimize…最后还是自己写了个模拟退火…反而可以了..无语

今天写两篇好了, 免得过几天又要整理,…而且今天感觉爬虫复习了挺多的…用R感觉比python好用多了啊…

关于豆瓣电影排名与时间的关系分析

从最简单的开始!

接着上次的GIS包继续说

不仅仅是geo数据 ggmap还可以绘制任何地方的地图

下面通过绘制近期中国地震地缘分布来练习一下使用
不多说上代码

library(ggmap)